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ABSTRACT 

A  boundary  integral  formulation  for  the  analysis  of 
circular  plate  bending  under  lateral  loads  is  developed  using 
Green's  functions.   The  formulation  specifically  applies  to 
annular  plates  with  arbitrary  boundary  conditions.   The  plate 
bending  solution  in  the  plastic  range  is  determined  using  a 
numerical  method  of  incremental  loading.   A  computer  program  to 
perform  the  required  calculations  was  developed  and  is 
presented.   Results  for  three  case  studies  are  included  and 
compared  with  results  obtained  by  other  methods.   Plate 
behavior  in  the  elastic  range  is  in  excellent  agreement  with 
other  analytical  solutions,  and  in  the  plastic  range  is  in 
reasonable  agreement  with  published  results  obtained  using  a 
finite  element  method. 
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NOMENCLATURE 


a  -  plate  outer  radius 

b  -  plate  inner  radius 

c  -  plate  outer  radius  for  Green's  function 

D  -  flexural  rigidity 

d  -  plate  inner  radius  for  Green's  function 

E  -  Young's  modulus 

h  -  plate  thickness 

I  -  area  moment  of  inertia  of  beam  cross  section 

Kr  -  radial  curvature 

Kg  -  tangential  curvature 

L  -  beam  length 

M,Mq  -  elastic  moments  for  beam  bending  problem 

Mr,M(jr  -  elastic  radial  moments  per  unit  length 

MrP  -  plastic  radial  moment  per  unit  length 

Mr  -  actual  radial  moment  per  unit  length 

Mq>Mg8  ~  elastic  tangential  moments  per  unit  length 

MqP  -  plastic  tangential  moment  per  unit  length 

Q,Qq  -  shears  for  beam  bending  problem 

Qr'^Gr  ~  elastic  radial  shears  per  unit  length 

Qr  -  actual  radial  shear  per  unit  length 

q,qQ  -  distributed  loads  per  unit  area 
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radial  radius  of  curvature 

tangential  radius  of  curvature 

radial  distance  from  plate  center 

radial  position  of  ring  load 

yield  stress 

deflections 

distance  along  beam 

lateral  distance  from  plate  center 

plastic  strain  increment 

total  radial  strain 

elastic  radial  strain 

plastic  radial  strain 

total  tangential  strain 

elastic  tangential  strain 

plastic  tangential  strain 

yield  strain 

slopes 

plate  boundary 

angular  position  in  tangential  direction 

biharmonic  operator 

loads  per  unit  length  for  beam  bending  problem 

magnitude  of  Green's  function  ring  load 

equivalent  stress 

radial  stress 

tangential  stress 

Poisson ' s  ratio 

plate  domain 
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CHAPTER   1 


INTRODUCTION 


The  solution  of  plate  bending  problems  in  the  plastic 
range  has  been  approached  using  finite  element,  and  more 
recently,  boundary  integral  methods  (also  called  boundary 
element  methods).   Unlike  the  finite  element  method,  which 
discretizes  the  inside  of  the  plate  into  a  number  of  small 
elements,  the  boundary  integral  method  discretizes  only  the 
plate  boundary,  thereby  reducing  the  order  of  the  problem  by 
one . 

Many  different  approaches  have  been  used  to  formulate 
plate  bending  problems  using  boundary  integrals,  with  perhaps 
the  most  attractive  proposed  by  Sternal'*  and  Bezine^>.   These 
authors,  working  independently,  developed  a  direct  boundary 
integal  method  using  Green's  functions  to  solve  the  general 
plate  bending  problem  in  the  elastic  range.   This  work  has  been 
extended  by  Moshaiov  and  Vorus^)  to  include  plastic  behavior. 

Symmetry  is  commonly  used  to  reduce  the  size  of  the  system 
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of  equations  needed  to  analyze  a  symmetric  structure.   In  the 
case  of  axisymmetr ical ly  loaded  circular  plates  the  problem 
becomes  one-dimensional.   For  the  linear  elastic  case,  it  has  a 
closed  form  analytic  solution.   However,  some  authors  have 
instead  treated  this  circular  plate  problem  as  a  two- 
dimensional  problem,  perhaps  for  the  sake  of  demonstration 
(Kamiya  and  Sawki^^).   They  have  used  symmetry  to  allow 
discretization  of  a  portion  of  the  boundary  instead  of  the 
entire  boundary. 

This  thesis  develops  a  general  formulation  for  solving 
annular  plate  bending  problems  using  boundary  integrals  that 
reduces  the  problem  to  essentially  a  one-dimensional  problem  by 
the  use  of  a  "ring"  type  Green's  function.   Such  a  formulation 
is  most  useful,  inspite  of  the  closed  form  solution,  as  it 
allows  the  analysis  of  plates  with  different  boundary 
conditions  and  loads  with  one  algorithm.   Moreover,  it  permits 
a  solution  in  the  plastic  range,  which  is  not  possible  with  a 
conventional  analytic  approach.   In  addition,  the  simplicity  of 
the  one-dimensional  formulation  given  here  has  a  unique 
significance  for  the  teaching  of  boundary  element  methods. 

Butterf ield<5>  has  developed  a  boundary  element 
formulation  for  the  one-dimensional  elastic  beam  bending 
problem.   Here,  a  similar  approach  is  taken  for  the  annular 
plate.   It  is  also  extended  to  include  non-linear  behavior 
using  an  incremental  load  method  as  outlined  by  Moshaiov  and 
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Vorus .   Numerical  results  for  three  different  annular  plate 
configurations  are  presented  and  compared  to  results  obtained 
using  other  methods.   A  discussion  of  these  results  follows,  as 
well  as  suggestions  for  future  work  in  this  area. 
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CHAPTER  2 


DEVELOPMENT  OF  THE  GOVERNING  DIFFERENTIAL  EQUATION 


This  chapter  reviews  the  derivation  of  the  governing 
differential  equation  for  a  circular  plate  subjected  to  a 
symmetrically  distributed  lateral  load.   The  derivation  is 
based  on  that  provided  by  Timoshenko<6> }  using  Kirchoff's 
assumptions.   It  is,  therefore,  applicable  only  to  the  linear- 
elastic  case.   Treatment  of  non-linear  behavior  is  addressed  in 
Chapter  5. 


2 . 1   Moment-Curvature  Relationships 

The  first  step  in  developing  the  governing  differential 
equation  for  a  circular  plate  is  to  find  a  relationship  between 
moment  and  curvature.   Curvature  in  the  radial  direction  (Kr) 
and  the  tangential  direction  (Kg)  for  a  symmetrically  loaded 
circular  plate  is  found  using  geometrical  considerations.   Kr , 
which  is  the  inverse  of  the  radius  of  curvature  in  the  radial 
direction  (Rr),  can  vary  as  a  function  of  the  distance  (r)  from 
the  plate  center.   Examining  the  plate  of  Figure  2-1  at  a  small 
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radial  distance  dr  away  from  r,  the  slope  (♦)  can  be  expressed 
in  terms  of  the  radial  distance  r  and  the  deflection  (w)  as 
follows 


dw 
dr 


(2-1-1) 


The  radial  radius  of  curvature  is 


dr 


d* 


giving  an  expression  for  the  radial  curvature  in  terms  of  w 


1 
R 


d* 
dr 


d  w 


dr 


(2-1-2) 
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Undeftec-ted 
position 


± 


Figure  2-1.   Deflection,  slope,  and  curvature  relationships 
in  circular  plate  bending 


Due  to  symmetry,  the  tangential  radius  of  curvature  is 
constant  at  any  given  radius  r,  and  for  small  deflections  can 
be  expressed  as 

r 

8    ♦ 

giving  an  expression  for  the  tangential  curvature 
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Q 


1  dw 
r  dr 


To  relate  moment  to  curvature,  consider  the  small  plate 
element  illustrated  in  Figure  2-2.   Define  Mr  and  Mg,  the 
radial  and  tangential  moments  per  unit  length  respectively,  as 
f ol lows : 


M 


M. 


h/2 
-h/2 
h/2 
-h/2 


a      z  dz 
r 


(2-1-3) 


a      z  dz 


where  h  is  the  plate  thickness.   Assuming  the  plate  is  thin  and 
hence  experiences  a  two-dimensional  stress  state,  Hooke's  law 
yields  the  following  expressions  relating  stresses  to  strains 


E 


(1-v2) 


9      f\7\ 
(  1-v  ) 


e         e 
r         0 


e         e 
e  _   +   v     e 
6  r 


(2-1-4) 


where  E  is  Young's  modulus,  v  is  Poisson's  ratio,  and  ere  and 
£ ge    are  the  elastic  radial  and  tangential  strains,  respec- 
tively.  In  Chapter  5,  which  discusses  non-linear  behavior,  a 
distinction  will  be  made  between  the  elastic  strain  and  the 
total  strain.   In  the  plastic  range,  the  total  strain  includes 
elastic  and  plastic  components,  in  accordance  with  the 
following  expressions: 
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e        p 
e    =     £     +   e 
r        r        r 


e   +   ,  P 
£9      ~         £9        +      £e 


where  the  superscript  p  refers  to  the  plastic  strain 
components.   Since  Chapter  2  treats  only  elastic  behavior,  the 
total  strain  is  simply  equal  to  the  elastic  strain. 


Figure  2-2.   Differential  circular  plate  element 


Substituting  Eq.  (2-1-4)  into  Eq.  (2-1-3)  gives 


19 


M 


(1-v  ) 


M 


(1-v  ) 


h/2 

[ 

-h/2 
ph/2 

[ 

-h/2 


£    +   V   £„      Z   dz 

r       fl 


(2-1-5) 


e _  +  v  £     z  dz 
8       r 


Finally,  from  geometrical  considerations,  strain  and  curvature 
are  related  by 


9 


K   z 
r 


K8  Z 


.  2 
d  W 


dr 


1  dw 


r  dr 


(2-1-6) 


which,  when  substituted  into  Eq.  (2-1-5),  give  the  desired 
moment-curvature  relationship  for  a  circular  plate 


M 


,  2 

d  w 


v  dw 
dr       r  dr 


(2-1-7) 


M 


-  D 


1  dw 


r  dr 


+   v 


,  2 

d  w 


dr 


where  D,  the  flexural  rigidity,  is  given  by 


E  h 


(2-1-8) 


12  (1-v  ) 
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2  .  2   Equilibrium  Equations 

For  the  moment-curvature  relationships  to  be  useful,  an 
equilibrium  equations  in  terms  of  moments  and  shears  must  be 
found.   The  governing  differential  equation  is  then  developed 
by  substituting  the  expressions  relating  moment  to  curvature 
into  the  equilibrium  equation. 

To  find  an  equilibrium  equation,  consider  first  the 
laterally  loaded  circular  plate  element  illustrated  in  Figure 
2-3,  where  Qr  is  the  radial  shearing  force  per  unit  length  and 
q  is  the  distributed  load  (force  per  unit  area).   The  couples 
acting  along  edges  ab  and  cd  due  to  radial  bending  moments  are 


ab: 


Mr  r  dS 


cd 


M   + 
r 


dM 


dr 


dr 


r  +  dr 


d8 


The  shear  forces  acting  along  the  edges  are 


ab 


Q  r  d8 


cd 


Q   + 
r 


dQ 


dr 


dr 


r  +  dr 


d9 


Q  + 

r 

i 

fcdr 

dr 

,,/' 

1-   1 

[■TT 

t 

lMr+^dr(  'W>. 

Wmw^ 

h 

*  Or 

t 

WI                                — 

•                               i                                * 

Z    1 

Figure  2-3 


Shear  forces  and  moments  acting  on  dif  f  erent  ia_i 
plate  element . 


Because  of  symmetry,  there  are  no  shear  forces  acting 
along  edges  ad  and  be,  leaving  only  the  tangential  bending 
moments 


ad: 
be: 


Ma   dr 
0 


These  couples  can  be  resolved  into  components  acting 
perpendicular  to  the  axis  ro  (which  are  equal  and  opposite  and 
thus  cancel),  and  components  acting  along  the  ro  axis  which 
have  a  total  magnitude  of 
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MQ  dr  d9 


Summing  up  all  the  couples  and  neglecting  the  small  change  in 
shear  force  across  the  element  (i.e.  dropping  the  dQ  terms) 
gives  the  following  equilibrium  equation 


Mr   + 


dM 


dr 


dr 


r  +  dr  J 


r  +  dr    dfl   -   M   r  d8 


-   Ms  dr  de   +   Q   r  dfl  dr    =    0 
o  r 


This  equation  is  further  simplified  by  ignoring  higher  order 
terms  to  give 


M    + 

r 


dM 


dr 


r       M„   +   Q   r 
6       r 


(2-2-1) 


Also,  from  equilibrium  in  the  z  direction 


dQ 


(2-2-2) 


dr 


2 . 3   Differential  Equation  for  Circular  Plate  Bending 

Having  found  the  desired  equilibrium  equations,  the  moment 
curvature  relationships  of  Section  2.1  are  now  used  to  develop 
the  governing  differential  equation  for  the  circular  plate 
bending  problem.   Specifically,  substituting  Eq.  (2-1-7)  into 
Eq.  (2-2-1)  results  in  a  third  order  differential  equation  to 
describe  the  response  of  a  circular  plate  to  lateral  loading 
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,  3 

d  w 


dr 


1  d  w 


r  dr 


1   dw 


r   dr 


Q 


(2-3-1) 


which  can  be  rewritten  to  provide  an  expression  for  shear 


Q 


,3 

d  w 


dr 


1  d  w 

+   —  

r  dr 


1   dw 
r   dr 


(2-3-2) 


Differentiating  Eq.  (2-3-1)  with  respect  to  r  and  using  Eq.  (2- 
2-2)  results  in  a  fourth  order  differential  equation  in  terms 
of  r,  w,  D,  and  q 


d  w 


dr 


2  d  w 

+   —  

,   3 

r  dr 


1   d  w 


1   dw 


r   dr 


+   — 


r   dr 


q 

D 


(2-3-3) 


which  can  be  expressed  in  a  more  concise  form  using  the 
biharmonic  operator  as 


4 

V  w 


q 

D 


(2-3-4) 


where 


dr 


2  d 

+   —  — 


1   d 


r  dr 


2   .   2 

r   ar 


1   d 

~3  ~~ 
r   dr 


(2-3-5) 
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CHAPTER  3 


A  BOUNDARY  INTEGRAL  METHOD 


This  chapter  describes  a  boundary  integral  method,  and 
how  it  can  be  used  to  provide  a  general  method  of  solution  of 
circular  plate  bending  problems.   A  one-dimensional  beam 
bending  problem  is  used  to  illustrate  this  method. 

3 . 1   General  Formulation 

Solution  of  Eq.  (2-3-3)  analytically  as  a  boundary  value 
problem  by  applying  known  boundary  conditions  is  possible 
( Timoshenko ) .   However,  this  approach  requires  special 
attention  for  each  plate  configuration  and  load  distribution. 
This  difficulty  is  compounded  when  considering  problems 
involving  plastic  behavior  which  should  be  solved  numerically 
due  to  the  non-linear  behavior  in  the  plastic  range. 

Therefore,  a  general  method  is  sought  wherein  the  same 
procedure  can  be  used  to  solve  a  variety  of  problems,  thus 
lending  the  problems  to  computer  solution  techniques.   One  such 
method,  adopted  for  this  work  to  examine  circular  plate 
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bending,  is  the  boundary  integral  method. 

Stern  and  Bezine,  working  independently,  obtained  a 
general  formulation  for  plate  bending  problems  in  terms  of 
boundary  integrals.   These  integrals  involve  displacements, 
slopes,  bending  moments,  and  shears  on  the  plate  boundary. 
Both  authors  used  the  following  reciprocal  identity,  which  can 
be  derived  using  Green's  identity: 


w.  7  w   -   w  7  w 


da 


QQ  w   +   MQ  ♦   - 


*G  M   -   wQ  Q   )  dr    (3-1-1) 


where 


a 

r 

D 
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w,wG   - 

*,*G   " 

M,MG   - 

Q,Qg  " 


plate  domain 

plate  boundary 

flexural  rigidity 

biharmonic  operator 

deflect  ions 

slopes 

bending  moments  per  unit  length 

shear  forces  per  unit  length 


To  arrive  at  the  boundary  integral  equation,  wG  is  selected  as 
a  Green's  function. 


3 . 2   Beam  Bending  Problem 

But t erf ield  has  presented  a  reduced  form  of  Eq.  (3-1-1) 
for  the  one-dimensional  case  of  beam  bending,  which  is 
developed  in  the  discussion  that  follows.   As  a  first  step, 
Butterfield  multiplies  both  sides  of  the  governing  differential 
equation  for  a  beam  by  a  second  function  and  integrates  over 
the  beam's  length,  giving 


EI  7  w  w_  dx 
u 


0 


\   w_  dx 


(3-2-1) 


0 


where  \    is  the  distributed  load  per  unit  length,  I  is  the  area 
moment  of  inertia  of  the  beam  cross  section,  and  the  operator 
V4  is  defined  for  the  beam  case  as 


dx 

Integrating  the  left  hand  side  of  Eq.  (3-2-1)  by  parts 
several  times,  Butterfield  develops  the  equation 


\    w_   -   w  EI  7 
G 


4wQ  )  dx 


-wQ  Q   +   #Q  N   -   MQ  ♦   +   QG  w 


Recognizing  that  for  the  second  function 


E  I  7  w. 


the  preceding  equation  can  be  rewritten  as 
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(  K    WG   "   W  *G  ) 


dx 


,L 


-wQ  Q   +   *Q  M   -   MQ  ♦   +   QQ  w 


(3-2-2) 


To  gain  some  physical  insight  into  Eq.  (3-2-2)  it  is 
useful  to  recall  Betti's  reciprocal  work  theorem.   This  theorem 
states  that  for  an  elastic  structure  subjected  to  two 
independent  causes  (e.g.  loads),  the  total  work  done  by  the 
first  cause  in  moving  through  the  displacements  resulting  from 
the  second  cause  is  equal  to  the  total  work  done  by  the  second 
cause  in  moving  through  the  displacements  produced  by  the  first 
cause^\   The  theorem,  as  it  applies  to  our  beam  bending 
problem,  can  be  better  understood  by  examining  the  examples  of 
Figure  3-1. 


28 


a, 


'999999/ 


Jl 
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♦ji 


Figure  3-1.   Examples  of  Betti's  reciprocal  work  theorem 

In  Figure  3-la,  the  product  of  load  P^  and  displacement 
5 i j  is  equal  to  the  product  of  load  Pj  and  displacement  8ji- 
Similarly,  in  Figure  3-lb,  the  product  of  moment  M^  and  slope 
♦ij  is  equal  to  the  product  of  moment  Mj  and  slope  *ji. 

This  concept  can  be  applied  to  understand  Eq.  (3-2-2), 
with  the  i  cause  the  actual  loading  or  applied  moment  condition 
for  the  system  of  interest  and  the  j  cause  the  loading  or 
applied  moment  condition  of  the  second  function.   In  this 
equation,  the  K    term  represents  the  external  load  in  the  system 
of  interest,  which  when  multiplied  by  the  deflection  wq 
described  by  the  second  function,  gives  a  reciprocal  work  terra. 
Similarly,  the  *q  term  is  a  slope  for  the  second  function, 
which  when  multiplied  by  the  moment  M  for  the  system  of 
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interest  yields  another  work  term.   This  concept  applies  to  all 
the  terms  of  the  equation,  which  are  combined  according  to 
Betti's  theorm  of  reciprocal  work  to  yield  the  given  equality. 

For  Eq.  (3-2-2)  to  be  useful  for  solving  the  beam  bending 
problem,  Butterfield  chooses  as  the  second  function  a  Green's 
function  which  is  based  on  a  singular  loading  condition;  that 
is  to  say,  a  point  load.   If  the  point  load  of  unit  magnitude 
is  applied  at  either  x=0  or  x=L,  the  second  intergral  term  on 
the  left  hand  side  of  Eq.  (3-2-2)  takes  on  the  value  of  the 
boundary  condition  of  the  beam  of  interest.   Usually,  some 
attention  should  be  given  to  cases  where  the  integrals  become 
s  ingular . 

This  yields  two  equations.   However,  in  general,  beam 
bending  problems  involve  four  unknown  boundary  conditions.   For 
example,  if  simple  supports  are  specified  on  both  ends,  the 
moments  and  deflections  at  the  two  ends  are  known  to  be  zero, 
while  the  shears  and  slopes  are  not  known.   Therefore,  two 
additional  equations  are  required. 

To  obtain  two  more  equations,  derivatives  of  the  first  two 
equations  are  taken.   This  gives  an  additional  Green's 
function.   Boundary  conditions  can  now  be  obtained  by  solving 
the  four  equations  simultaneously.   Knowing  the  boundary 
conditions,  the  beam  problem  can  be  solved  for  interior  points 
using  Eq.  (3-2-2)  and  appropriate  derivatives  with  the  point 
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load  positioned  at  the  location  of  interest. 


The  formulation  for  the  beam  bending  problem  outlined 
above  is  a  direct  approach  based  on  Green's  identity. 
Butterfield  also  develops  a  boundary  integral  equation  using 
the  more  intuitive  indirect  method.   For  the  circular  plate 
formulation  that  follows  in  Chapter  4,  only  the  direct  method 
is  used. 


CHAPTER  4 


DEVELOPMENT  OF  THE  FORMULATION  FOR  CIRCULAR  PLATES 


This  chapter  presents  a  general  method  for  solving 
circular  plate  bending  problems  in  the  elastic  range  using 
boundary  integrals.   The  resulting  integral  equations  are 
analagous  to  those  developed  by  Butterfield  for  the  one- 
dimensional  beam  problem,  which  are  discussed  in  Chapter  3. 
Treatment  of  non-linear  behavior  is  discussed  in  Chapter  5. 

4 . 1   Description  of  Plate  Geometry 

Figure  4-1  depicts  a  symmetrically  loaded  annular  plate, 
and  includes  notation  that  will  be  referred  to  throughout  the 
development  of  the  integral  equation  and  the  discussion  of  the 
selected  Green's  functions. 

Important  symbols  are  defined  as  follows: 


w,wG 

#,*G 

Mr>M(jr      - 

M8»MG8 

Qr.QGr      - 

Qq,QG8      - 
a 

b 

h 

r 

q 


deflect  ions 

slopes 

radial  bending  moments  per  unit  length 

tangential  bending  moments  per  unit  length 

radial  shear  forces  per  unit  length 

tangential  shear  forces  per  unit  length 

plate  outer  radius 

plate  inner  radius 

plate  thickness 

radius  of  interest 

distributed  load  per  unit  area 
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\^S^ 


Mo- 


a 
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Figu 


Symmetrically  loaded  annular 
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4 . 2   Development  of  Integral  Equation 

Though  the  laterally  loaded  circular  plate  bending  problem 
is  two-dimensional,  symmetry  reduces  the  problem  to  essentially 
a  one-dimensional  type  problem.   Therefore,  it  is  expected  that 
the  integral  equation  developed  by  Butterfield  for  the  beam 
case  can  be  adapted  to  the  case  of  the  circular  plate.   This 
equation  (Eq.  (3-2-2))  is  repeated  below: 


(  K    WG   "   W  \3  )  dx 


-wQ  Q   +   *Q  M   -   MQ  i      +   QQ  w 


0 


To  obtain  some  indication  of  what  the  circular  plate 
integral  equation  will  look  like,  a  non-rigorous  method  is 
first  used  to  adapt  Eq.  (3-2-2)  to  the  circular  plate  case, 
followed  by  a  more  rigorous  mathematical  development. 

In  the  non-rigorous  approach,  a  first  step  is  to  replace 
the  integral  on  the  left  hand  side  with  an  area  integral  and 
the  line  load  per  unit  length  \    with  the  load  per  unit  area  q. 
Next,  adjustments  are  made  to  reciprocal  work  terms  on  the 
boundaries  (right  hand  side  of  Eq.  (3-2-2)).   Specifically,  the 
shear  terms  are  replaced  with  the  product  of  the  radial  shear 
per  unit  length  and  the  total  length  (2ur) ,  observing  that  for 
symmetrical  loading,  the  tangential  shear  Qg  is  0.   Moment 
terms  are  replaced  with  radial  moments,  since  there  are  no 
tangential  moments  on  the  boundaries.   These  moments  must  also 
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n 


be  multiplied  by  the  length  (2ir)  to  give  total  moment.   The 
resulting  expression  is 


,2i 


J0   J 


q  w    -   w  q     r  dr  dfl 


2"r  (  "WG  Qr  +   »G  Mr  "   MGr  *  +      QGr  W  ) 


(4-2-1) 


which  reduces  to 


(  q  wQ   -   w  qQ  ) 


r  dr 


"WG  Qr  +   *G  Mr  "   MGr*   +   QGrW  ) 


(4-2-2) 


So  far,  it  has  only  been  asserted  that  the  integral 
equation  for  the  circular  plate  case  should  have  a  form  similar 
to  Eq.  (4-2-2).   To  obtain  the  exact  integral  equation,  a 
direct  method  is  used. 

As  a  starting  point  for  the  direct  method,  consider  the 
following  expression: 


U 


a  A2  A2 

d  w  d  w 
(j 

'a        *  A        2 

,  dr  dr 


r  dr 


This  expression  can  be  written  slightly  differently  as  V 


35 


3  A2  A7 

d  w   d  w, 


2        2 

dr    dr 


r  dr 


Integrating  both  expressions  by  parts  yields 


dw_  d  w 

U    =    — r 

dr   dr 


dw„  [    1   d  w 


dr 


d  w  \ 


r   dr 


dr 


r  dr 


and 


dw  d  w 
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dr  dr 


dw 
dr 


1   d  w. 


r   dr 


,3 

d  w. 


dr 


r  dr 


Recalling  Eq.  (2-3-2),  it  is  possible  to  rewrite  these 
expressions  in  terms  of  shear  Qr 


U 


dw„  d  w 
dr   dr 


,a 


dw   Q 

G   r    , 
—  —   r  dr 

dr   D 


1  dw_  dw 
u 

r  dr   dr 


dr 
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dw  d  w, 


dr  dr 
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dw  QQr 

—  —    r  dr 
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1  dw  dw. 


r  dr  dr 


dr 


Setting  U  equal  to  V  and  canceling  some  terms  gives 
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dw   d  w 
G 

dr   dr 


dw   Q 

G   r    , 
—  —   r  dr 

dr   D 


dw  d  w, 


dr  dr 


Jb 


dw   Q 

—  —    r  dr 

dr   D 


The  two  integral  terms  above  are  integrated  by  parts,  with 
the  operator  of  Eq.  (2-3-5)  used  to  simplify  the  expressions 


dwQ  Qr 

—   —   r  dr 

dr   D 
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G    D 
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which,  when  substituted  into  the  previous  equation,  yield 
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After  some  simplification,  this  expression  becomes 
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The  last  equation  bears  a  close  resemblance  to  Eq.  (4-2- 
2),  except  for  sign  differences  and  the  fact  that  terms 
involving  the  second  derivative  of  deflection  have  yet  to  be 
replaced  by  terms  involving  moment.   To  accomplish  the  latter, 
the  following  expression  is  added  to  the  third  term  of  the 
right  hand  side  of  the  preceding  expression,  and  subtracted 
from  the  fourth  term 


v  dw  dw 
r  dr  dr 


G 


Through  Eq.  (2-1-7)  this  gives 


V  w  r  w,,   -   7 


4wQ  r  w  ) 


dr 


WG  Qr   "   W  QGr  +   *  MGr 


♦  _  M 
G   r 


(4-2-3) 


which  can  be  rewritten  using  Eq.  (2-3-4)  as 


Jb 


q  wQ   -   qQ  w    r  dr 


WG  Qr   "   W  QGr  +   *  MGr  "   *G  Mr 


(4-2-4) 


Remarkably,  with  the  exception  of  a  sign  difference  that  is 

attributable  to  the  difference  in  sign  convention,  this 

expression  is  identical  to  Eq.  (4-2-2),  which  was  asserted 

based  on  the  integral  equation  developed  by  Butterfield. 
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4. 3   Selection  of  a  Green's  Function 

As  is  the  case  for  the  beam  described  in  Chapter  3,  Eq. 
(4-2-4)  is  useful  only  if  wg  is  chosen  as  a  Green's  function 
with  certain  properties.   As  will  become  evident  later  in  this 
chapter,  the  Green's  function  must  describe  a  ring  loaded  plate 
to  solve  circular  plate  bending  problems,  just  as  a  Green's 
function  describing  a  point  loaded  beam  was  used  to  solve  the 
beam  bending  problem. 

To  solve  the  beam  problem,  Butterfield  choses  a  Green's 
function  describing  an  infinitely  long  beam  with  a  point  load. 
It  is  not  essential  that  a  Green's  function  for  an  infinite 
beam  be  used.   In  fact,  as  one  would  expect  from  Betti's 
reciprocal  work  theorem,  a  Green's  function  describing  a  finite 
geometry  will  work  as  well,  provided  the  function  also 
describes  the  response  to  a  point  load. 

Applying  But terf ield ' s  method  to  the  circular  plate  case, 
a  Green's  function  representing  the  ring  loaded  plate  shown  in 
Figure  4-2  is  adopted.   The  Green's  function  and  associated 
derivatives  are  included  in  Appendix  A.   It  should  be  noted 
that  the  outer  and  inner  radii  of  the  Green's  function  (c  and 
d)  need  not  coincide  with  the  outer  and  inner  radii  of  the 
actual  plate  (a  and  b). 
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Figure  4-2.   Plate  geometry  for  the  selected  Green's 
f unct  ion 


That  a  simple  support  is  chosen  for  the  outer  plate  radius  is 
not  important;  a  plate  clamped  at  both  the  inner  and  outer 
radii  would  also  have  been  suitable.   What  is  important  is  that 
the  loading  condition  for  the  Green's  function  case  be  a  ring 
load.   To  simplify  calculations,  a  ring  load  of  magnitude  1  is 
used. 

4 . 4   Solution  of  Integral  Equation 

Having  selected  a  Green's  function,  the  next  step  is  to 
obtain  a  form  of  Eq.  (4-2-4)  suitable  to  solve  the  actual  plate 
bending  problem.   Rearranging  terms,  Eq.  (4-2-4)  can  be 
rewritten  as 
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q„  w  r  dr 
u 


-WG  Qr  +   *o  Mr  -   MGr  *   +   QGr  w 


q  w   r  dr 


(4-4-1) 


The  ring  load  may  be  expressed  in  terms  of  the  dirac  delta 
function  8 ( r , r0) 


D  V  w, 


5(r,ro) 


Substituting  this  expression  into  Eq.  (4-4-1)  and  making  use  of 
the  following  property  of  a  dirac  delta  function 


8(r,r  )  w  r  dr 
o 


r   w(r=r  ) 
o       o 


gives 


r   w( r=r  ) 
o       o 


-WG  Qr  +   *G  Mr   "   MGr  *      +   V  " 


q  w_  r  dr 


(4-4-2) 


By  choosing  r0  at  radius  a,  and  then  at  radius  b,  the 
deflection  w  in  the  left  hand  side  of  Eq.  (4-4-2)  is  the 
deflection  at  the  boundary,  which  is  either  a  known  or  unknown 
boundary  condition.   Hence,  two  equations  are  available  to 
solve  for  the  four  unknown  boundary  conditions. 
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To  obtain  two  more  equations,  Eq.  (4-4-2)  is 
differentiated  with  respect  to  r0.   This  yields  the  second 
Green's  function,  which  is  evaluated  with  r0  at  a  and  then 
again  at  b.   The  second  Green's  function  corresponds  to  the 
slope  of  the  plate,  which  at  r=a  and  r=b  is  again  either  a 
known  or  an  unknown  boundary  condition.   There  is  thus  now  a 
complete  set  of  four  equations  to  solve  for  the  four  unknown 
boundary  conditions. 

Eq.  (4-4-2)  can  be  rewritten  as 


w( r=r  ) 
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"WG  Qr  +   *G  Mr   "   MGr  *   +   QGr  K 


+       — 


q  w„  r  dr 


(4-4-3) 


Differentiating  this  expression  with  respect  to  r0  gives 
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(4-4-4) 


Eqs.  (4-4-3)  and  (4-4-4),  evaluated  both  for  r0=a  and  r0=b, 
give  rise  to  a  system  of  four  equations  that  when  solved 
simultaneously,  yield  the  four  unknown  boundary  conditions, 
matrix  form,  these  equations  can  be  written  as 


In 
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(4-4-5) 


Using  Gauss  elimination,  the  preceding  system  of  equations  is 
solved  for  the  four  unknown  boundary  conditions. 

The  final  step  is  to  calculate  deflection,  slope,  moment 
and  shear  across  the  plate.   Deflection  and  slope  are  found 
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using  Eqs.  (4-4-3)  and  (4-4-4)  by  substituting  in  the 
calculated  boundary  conditions  and  evaluating  the  expressions 
at  the  value  of  r0  of  interest.   Moment  and  shear  are 
determined  using  Eqs.  (2-1-7)  and  (2-3-2),  which  require 
calculation  of  the  second  and  third  derivatives  or  w  with 
respect  to  r0.   Differentiating  Eq.  (4-4-4)  results  in  the 
desired  expressions 


d  w(r=r  ) 

7? 


.  2 

d  w 
dr2    r 


dr 


d2M 


Gr 


d2Q 


*   + 


Gr 


w 


dr 


dr 


i  a 


+    — 
r 


a    a2 
d  w. 


r  dr 


dr 
o   b      o 


2   dw(r=r  ) 


r   dr 
o    o 


(4-4-6) 


and 
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making  a  complete  solution  for  the  circular  plate  bending 
problem  in  the  elastic  range  possible. 
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CHAPTER  5 


NON-LINEAR  BEHAVIOR 


This  chapter  develops  a  boundary  integral  formulation  for 
circular  plate  bending  problems  that  includes  non-linear 
material  behavior.   Specifically,  the  boundary  integral 
equation  developed  in  Chapter  4  for  the  elastic  plate  is 
modified  to  include  plastic  moment  terms  to  account  for  the 
plastic  behavior.   An  incremental  load  method  similar  to  that 
used  by  Moshaiov  and  Vorus  can  then  be  applied  to  arrive  at  a 
solution  to  the  plate  bending  problem  once  yielding  has  begun. 

5 . 1   Plastic  Stress-Strain  Relationships 

When  considering  plasticity,  strain  can  be  thought  of  as 
having  two  components:  an  elastic  component  and  a  plastic 
component.   In  cylindrical  coordinates,  which  is  applicable  to 
the  two-dimensional  circular  plate  case,  the  total  strain  e    can 
be  expressed  in  terms  of  these  two  components  as 


e       p 

£      =      £       +    t 

r        r        r 


e        p 


(5-1-1) 


where  the  superscripts  e  and  p  refer  to  the  elastic  and  plastic 
strain  components,  respectively. 

Eq.  (2-1-4),  which  provides  the  relationship  between 
stress  and  strain  for  a  linearly  elastic  material,  can  thus  be 
rewritten  in  terms  of  the  total  strain  and  plastic  strain  as 
f ol lows : 
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5 . 2   Plastic  Moment-Curvature  Relationships 

To  develop  the  moment-curvature  relationships  with 
plasticity,  Eq.  (5-1-2)  is  substituted  into  Eq.  (2-1-3)  to  give 
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These  expressions  can  be  rewritten  a! 
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where  the  radial  and  tangential  plastic  moments,  MrP  and  MgP, 
are  defined  as 
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Recalling  the  relationships  for  total  strain  of  Eq.  (2-1-6)  and 
performing  the  integration  over  the  plate  thickness  yield  the 
moment-curvature  relationship  for  a  circular  plate  with 
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5 . 3   Governing  Differential  Equation  With  Plasticity 

Chapter  2  develops  the  governing  differential  equation  for 
a  circular  plate  in  the  elastic  range  by  substituting  the 
elastic  moment-curvature  relationships  into  the  elastic 
equilibrium  equation.   In  this  section,  a  similar  approach  is 
used  to  develop  the  governing  differential  equation  with 
plasticity.   The  difference  here  is  that  the  moment-curvature 
relationships  (Eq.  (5-2-3))  include  plastic  terms. 

Substituting  Eq.  (5-2-3)  into  the  equilibrium  equation  for 
a  circular  plate  (Eq.  (2-3-1))  gives 
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where  the  tilda  is  used  to  distinguish  the  fact  that  the  shear 
includes  plastic  effects.   Differentiating  this  expression  once 
with  respect  to  r  gives  the  fourth  order  governing  differential 
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equation  for  the  plate  with  plasticity 
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which  can  be  rewritten  using  the  operator  defined  in  Eq.  (2-3 
5 )  as 
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This  equation  is  identical  to  the  equilibrium  equation  for  the 
elastic  case  with  the  exception  of  the  terms  that  include 
derivatives  of  the  plastic  moments.   Note  that  these  terms  may 
be  thought  of  as  additional  pseudo  loads  that  must  be  combined 
with  the  actual  load  to  account  for  the  plastic  behavior. 

5 . 4   Boundary  Integral  Formulation  With  Plasticity 

In  the  preceding  chapter,  the  governing  differential 
equation  is  used  to  develop  a  boundary  integral  formulation  for 
a  circular  plate  in  the  elastic  case.   This  section  presents  a 
similar  approach  to  develop  a  boundary  integral  formulation 
for  the  plastic  case. 


As  a  starting  point,  Eqs.  (2-3-4)  and  (5-3-2)  are 
substituted  into  Eq.  (4-2-3)  to  give 
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This  equation  is  not  yet  in  a  form  that  is  usable  for 
solving  the  circular  plate  bending  problem.   In  the  first 
place,  there  are  several  terms  that  include  derivatives  of  the 
radial  or  tangential  plastic  moments,  which  are  not  generally 
known  across  the  plate.   Secondly,  the  shear  and  moment  terms 
Qr  and  Mr  evaluated  at  the  boundaries  are  not  the  actual  shears 
and  moments,  since  they  were  developed  based  on  elastic 
considerations  only. 

To  resolve  these  concerns,  the  integral  term  involving 
plastic  moments  in  Eq.  (5-4-1)  is  integrated  by  parts  as 
f ol lows : 
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Integrating  by  parts  a  second  time  and  collecting  terms  gives 
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Using  the  above  relationship  and  choosing  a  Green's  function 
that  describes  a  ring  loaded  plate  as  is  done  for  Eq.  (4-4-2), 
Eq.  (5-4-1)  becomes 
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Recall  that  Qr  and  Mr  represent  shear  and  moment  for  the 
elastic  case.   However,  they  are  not  the  actual  shear  and 
moment  for  the  plastic  case.   With  plasticity,  the  actual 
shears  and  moments  are  combinations  of  elastic  and  plastic 
terms.   From  Eq.  (5-3-1),  the  shear  relationship  with 
plasticity  is 
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and  from  Eq.  (5-2-3),  the  moment  relationship  with  plasticity 


is 
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where,  as  with  the  case  for  shear,  the  tilda  is  used  to 
distinguish  the  fact  that  the  moment  term  includes  plastic 
effects.   When  substituted  into  the  preceding  expression,  these 
expressions  give  the  desired  boundary  integral  formulation  for 
the  circular  plate  bending  problem 
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By  evaluating  Eq.  (5-4-2)  with  the  ring  load  positioned 
both  at  r0=a  and  r0=b,  two  equations  are  obtained  to  solve  for 
the  unknown  boundary  conditions.   Eq.  (5-4-2)  is  differentiated 
with  respect  to  r0  to  obtain  a  second  Green's  function.   This, 
via  Eq.  (2-1-1),  gives  an  expression  for  slope  of  the  actual 
plate.   The  second  Green's  function  is  evaluated  with  the  ring 
load  positioned  at  r0=a  and  r0=b,  to  provide  the  remaining  two 
equat  ions . 

Differentiating  Eq.  (5-4-2)  with  respect  to  r0  gives 
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In  matrix  form,  the  system  of  equations  to  be  solved  for 
the  four  unknown  boundary  conditions  for  the  non-linear  case 
can  be  written  using  the  notation  of  Eq.  (4-4-5)  as  follows: 


54 


aw(  a) 
bw(b) 
■a#(a) 
-b#(b) 


=  a 


-wQ(a,a)  *G(a,a)   -MQr(a,a)   QQr(a,a) 

-wQ(a,b)  #Q(a,b)   -MQr(a,b)   QQr(a,b) 

-wG(a,a)  *G(a,a) 

-w_(a,b)  *_(a,b) 


G 


-MQr(a,a)   QQr(a,a) 

-M*  (a,b)   q'  (a,b) 
Gr         Gr 


-wQ(b,a) 
-wQ(b,b) 
-wQ(b,a) 
-wG(b,b) 


#G(b,a)   -MQr(b,a) 


QGr(b,a) 


*G(b,b)   -MQr(b,b)   QQr(b,b) 


#G(b,a)   -MQr(b,a) 
#G(b,b)   -MQr(b,b) 


QGr(b,a) 
QGr(b,b) 


Qr(b) 

Mr(b) 

♦  (b) 

w(b) 


q  wQ(r,a)  r  dr 

ra 

q  wQ(r,b)  r  dr 

a 

q  Wg(r, a)  r  dr 

a 

q  wG(r,b)  r  dr 


dw  (r,a) 

—       % 
dr        9 


w(a) 
w(b) 


d  w  (r,a) 

r  M 

2  r 


dw_(r,b) 
dr 


dwQ(r,a)^ 


dr 


d  wQ(r,b) 


dr 


r  M 


dr 


dr 


M 


d  wQ(r,a) 
dT7 


Q  (a) 
r 

Mr(a) 
♦  (a) 
w(a) 


dr 


dr 


dwG(r,b)l 
dr 


d  wG(r,b)l 


M 


dr 


dr 


dr 


(5-4-4) 


r  M 


dr 


r  M 


dr 


55 


Observe  that  Eq.  (5-4-4)  is  identical  to  Eq.  (4-4-5) 
except  for  the  additional  moment  terms  in  the  last  matrix  of 
Eq.  (5-4-4).   These  additional  terms  account  for  the  plastic 
behavior . 

To  find  moments  and  shears  across  the  plate,  expressions 
for  the  second  and  third  derivatives  of  deflection  with  respect 
to  r0  are  required.   These  expressions  follow: 
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CHAPTER  6 


NUMERICAL  SOLUTION 


This  chapter  describes  the  numerical  procedure  that  was 
used  to  solve  the  circular  plate  bending  problem  in  both  the 
elastic  and  plastic  ranges.   The  incremental  load  approach  used 
in  the  plastic  range  is  first  described.   A  brief  description 
of  the  computer  program  that  was  developed  based  on  the 
equations  of  Chapters  4  and  5  follows,  including  a  discussion 
of  some  aspects  of  numerical  methods  used  in  this  program. 

6 . 1   Incremental  Load  Method 

Once  plastic  deformation  starts  to  occur  in  the  plate,  the 
linear  relation  between  stress  and  strain  is  not  valid.   As 
discussed  by  Mendelson^^,  strains  in  the  plastic  range  are  no 
longer  uniquely  determined  by  the  stresses,  and  in  fact  depend 
on  the  loading  history.   Therefore,  it  is  not  possible  to  use 
the  relationships  developed  in  Chapter  5  to  directly  arrive  at 
a  plate  bending  solution  given  a  loading  condition  outside  the 
elastic  range.   Rather,  an  incremental  approach  must  be  used 
wherein  the  load  is  increased  in  small  steps  once  yielding 
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occurs  at  any  point  in  the  plate,  and  the  complete  stress  state 
for  the  plate  is  determined  before  load  again  is  increased. 

The  incremental  load  method  used  for  this  analysis  is 
based  on  that  adopted  by  Moshaiov  and  Vorus,  and  is  frequently 
seen  in  finite  element  solutions  of  structural  analysis 
problems.   The  principal  steps  of  this  method  are  outlined 
below: 

1.  Using  the  plate  parameters  and  loading  conditions  for 
the  bending  problem  at  hand,  an  elastic  solution  is  obtained 
with  the  equations  developed  in  Chapter  4. 

2.  The  load  at  which  the  onset  of  yielding  will  occur  in 
the  plate  is  calculated  based  on  the  stress  state  determined  in 
Step  1.   The  selected  yield  criterion  is  the  von  Mises 
criterion,  which  for  the  symmetrically  loaded  circular  plate 
case  has  the  form: 


-  a       a  . 
r   9 


+  a 


where  ae    is  the  equivalent,  or  effective,  stress  which 
represents  the  von  Mises  yield  surface.   Yielding  will  occur 
when  <7e  equals  or  exceeds  the  uniaxial  yield  stress  (Sy). 

3.    Unknown  boundary  conditions  at  yielding  are 
determined  using  Eq.  (5-4-4),  with  plastic  moments  initially 


59 


set  equal  to  zero.   Total  strains  at  predetermined  integration 
points  are  calculated  using  Eqs.  (2-1-6),  (5-4-3)  and  (5-4-5). 
From  Eq.  (2-1-4),  the  stress  state  at  the  integration  points  at 
yielding  is  determined,  and  is  stored  to  be  later  updated  as 
load  is  increased.   Also  stored  for  later  updating  are  the 
deflections,  slopes,  moments  and  shears  across  the  plate. 

4.  The  load  which  caused  yield  is  then  increased  by  a 
small  incremental  amount.   The  increments  of  the  total  strains 
resulting  from  this  incremental  load  are  calculated  using  Eq. 
(2-1-6)  in  an  incremental  form,  from  which  elastic  strains  can 
be  calculated  by  means  of  Eq.  (5-1-1). 

5.  Using  Eq.  (2-1-4),  the  elastic  stress  increment 
corresponding  to  the  applied  incremental  load  is  calculated. 
This  incremental  stress  is  added  to  the  stress  stored  in  Step  3 
to  give  the  total  stress  at  the  end  of  the  load  increment.   A 
check  is  made  throughout  the  plate  using  the  yield  criterion  of 
Step  2  to  determine  where  yielding  has  occurred. 

6.  In  those  plate  regions  where  yielding  has  occurred, 
the  plastic  strain  increment  (that  is,  the  plastic  strain  due 
to  the  incremental  load)  is  calculated.   Only  materials  that 
exhibit  elas t ic-perf ect ly  plastic  behavior  as  depicted  in 
Figure  6-1  have  been  examined.   However,  the  method  can  be  used 
with  other  material  behavior,  such  as  strain  hardening.   For 
materials  exhibiting  elas tic-perfect ly  plastic  behavior,  the 
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plastic  strain  increment  8ep  is  the  total  incremental  strain 
determined  from  Eq.  (2-1-6).   This  is  shown  in  Figure  6-1. 

7.  The  plastic  moments  are  determined  using  Eq.  (5-2-2). 

8.  Steps  3  through  6  are  repeated,  with  the  plastic 
moments  calculated  in  Step  7  used  as  an  input  to  Eq.  (5-4-4). 
Step  7  is  repeated  and  the  plastic  moments  compared  to  the 
moments  calculated  previously  in  Step  7.   Iterations  are 
performed  as  necessary  until  these  values  converge,  according 
to  a  predetermined  convergence  criteria. 

9.  After  convergence  of  plastic  moments  is  achieved  in 
Step  8,  stresses,  plastic  strains,  deflections,  slopes,  moments 
and  shears  throughout  the  plate  are  updated. 


10.   Another  increment  of  load  is  applied,  and  Steps  3 
through  9  above  are  repeated,  using  the  plastic  moment  from  the 
previous  load  step  as  an  initial  estimate  for  plastic  moment. 
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Figure    6-1 


Uniaxial  stress-strain  curve  for  elastic- 
perfectly  plas-tic  material. 


6. 2   Computer  Program  Description 

Computer  programs  written  in  the  Fortran  77  language  have 
been  developed  to  solve  the  circular  plate  bending  problem  and 
are  included  as  Appendix  B.  Three  programs,  supported  by  nine 
subroutines,  are  used. 

The  programs  INPUT  and  LOAD  create  data  files  which  are 
used  by  the  program  MAIN  to  solve  the  bending  problem.   INPUT 
and  LOAD  prompt  the  user  for  information  describing  the  problem 
to  be  solved,  such  as  material  properties,  plate  geometry, 
boundary  conditions  and  plate  loading  conditions.   In  addition, 
certain  adjustable  parameters  are  input,  which  include  step 
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sizes  for  the  load  increments,  number  of  integration 
increments,  and  the  desired  percentage  for  convergence  of  the 
plastic  moment  increments.   The  program  MAIN  determines  the 
plate  bending  solution,  both  in  the  elastic  and  plastic  ranges 


A  flow  diagram  outlining  the  major  elements  and  overall 
logic  path  of  this  program  is  shown  in  Figure  6-2. 
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Figure  6-2.   Computer  program  flow  diagram 
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6 . 3   Numerical  Methods 

Eqs.  (5-2-1),  (5-4-4),  (5-4-5)  and  (5-4-6)  require 
integrations  to  be  performed  across  the  plate.   These 
integrations  are  accomplished  by  a  trapozoidal  rule  numerical 
integration  scheme.   It  was  found  that  10  increments  across 
both  the  plate  radius  and  the  plate  half  thickness  gave 
satisfactory  results. 

To  ensure  convergence  of  plastic  moments,  a  root  mean 
square  average  of  the  the  plastic  moments  at  each  station 
across  the  plate  is  calculated  and  compared  to  the  root  mean 
square  average  from  the  previous  iteration.   An  agreement  of 
less  than  0.1%  was  used  for  the  analysis  work  for  which  results 
appear  in  Chapter  7. 
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CHAPTER  7 


RESULTS  AND  DISCUSSION 


This  chapter  presents  results  of  analyses  performed  on 
three  different  plate  configurations  using  the  computer  program 
described  in  the  preceding  chapter.   Results  in  the  elastic 
range  for  all  three  cases  are  compared  with  the  analytic 
solution  from  Roark<^>)  an(j  for  one  case  in  the  plastic  range 
where  published  results  using  another  solution  method  are 
available.   Descriptions  of  the  case  studies  are  given  in 
Sections  7.1  through  7.3.   A  discussion  of  the  results  is  given 
in  Section  7.4. 

7 . 1   Simply  Supported  Elasto-Plast ic  Annular  Plate 

The  annular  plate  shown  in  Figure  7-1  is  first  examined. 
The  plate  is  simply  supported  at  both  the  inner  and  outer 
radii,  and  is  subjected  to  a  uniform  lateral  load.   A  material 
exhibiting  elast ic-perf ect ly  plastic  behavior  as  depicted  in 
Figure  6-1  is  used.   The  yield  stress  (Sy)  is  16  ksi  and  the 
Young's  modulus  is  lOxlO3  ksi.   Material  properties  and  plate 
dimensions  were  selected  to  allow  comparison  with  results 
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obtained  using  a  finite  element  method  by  Armen  et  al<^0>. 

Plate  deflection  at  yielding  is  shown  in  Figure  7-2,  which 
also  includes  results  from  the  analytic  solution  given  by  Roark 
for  comparison.   Plate  deflections  for  the  present  solution  at 
several  loads  in  the  plastic  range  are  shown  in  Figure  7-3,  as 
well  as  results  obtained  by  Armen  et  al  using  a  finite  element 
method  of  analysis. 

Figure  7-4  shows  the  plastic  zones  in  the  plate  at  the 
three  load  conditions  depicted  in  Figure  7-3.   Results  obtained 
by  Armen  et  al  are  also  included. 
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Figure  7-1.   Simply  supported  annular  plate  problem 
description . 
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ure  7-2.   Deflection  of  simply  supported  annular  pia    at 
the  onset  of  yielding. 
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Figure  7-3.   Deflection  of  simply  supported  annular  pial.  in 
the  plastic  range. 
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'igure  7-4.   Plastic  zones  in  simply  supported  annular  pit-.- 

7 . 2   Clamped  Elasto-Plast ic  Annular  Plate 

The  annular  plate  shown  in  Figure  7-5  is  next  examined. 
The  plate  is  clamped  at  both  the  inner  and  outer  radii,  and  is 
again  subjected  to  a  uniform  lateral  load.   A  material 
exhibiting  the  same  properties  as  described  in  Section  7.1  is 
used.   For  the  clamped  annular  plate  case,  results  obtained 
using  other  methods  could  not  be  found  for  the  plastic  range. 
Consequently,  material  properties  and  plate  dimensions  were 
selected  to  allow  a  qualitative  comparison  with  results  for  the 
simply  supported  plate  of  Section  7.1 

Plate  deflections  at  the  onset  of  yielding  are  shown  in 
Figure  7-6,  and  at  several  loads  in  the  plastic  range  in  Figure 
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7-7.   The  deflection  at  yielding  using  the  analytic  solution 
given  by  Roark  is  shown  in  Figure  7-6  for  comparison.   Finally, 
Figure  7-8  shows  the  plastic  zones  in  the  plate  at  the  three 
load  conditions  depicted  in  Figure  7-7. 
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Figure  7-5.   Clamped  annular  plate  problem  descriptic 
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Figure    7-6.       Deflection    of    clamped    annular    plate    at    the 
onset    of   yielding 
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Figure  7-7.   Deflection  of  clauped  annular  plate  in  the 
plastic  range 
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Figure  7-8.   Plastic  zones  in  clamped  annular  plate. 

7 . 3   Simply  Supported/Guided  Elas to-Plastic  Annular  Plate 

The  third  annular  plate  to  be  examined  is  shown  in  Figure 
7-8.   The  plate  is  simply  supported  at  the  outer  radius,  with 
the  inner  radius  guided.   The  plate  is  again  subject  to  a 
uniform  lateral  load.   A  material  exhibiting  the  same 
properties  described  in  Section  7-1  is  assumed. 

For  this  case,  results  obtained  using  other  methods  could 

again  not  be  found  for  the  plastic  range.  However,  a  solution 

for  a  continuous  simply  supported  circular  plate  has  been 

obtained  by  Moshaiov  and  Vorus,  and  offers  a  qualitative 

comparison  with  the  results  for  the  simply  supported/guided 
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annular  plate. 

Plate  deflection  at  the  onset  of  yielding  is  shown  in 
Figure  7-10  and  at  several  loads  in  the  plastic  range  in  Figure 
7-11  for  the  present  analysis  method.   The  deflection  at 
yielding  obtained  using  the  analytic  solution  from  Roark  is 
shown  in  Figure  7-10  for  comparison.   Figure  7-12  shows  the 
plastic  zones  in  the  plate  at  the  three  load  conditions 
depicted  in  Figure  7-11  for  the  present  solution.   Corres- 
ponding results  from  Moshaiov  and  Vorus  for  the  continuous 
plate  of  Figure  7-13  are  shown  in  Figures  7-14  and  7-15. 
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Figure  7-9.   Simply  supported/guided  annular  plate  probl  .3 
description . 
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Figure  7-11.   Deflection  of  simply  supported/guided  annular 
plate  in  the  plastic  range 
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7 . 4   Discussion  of  Results 

The  results  presented  in  Sections  7.1  through  7.3  for  the 
elastic  range  are  uniformly  in  excellent  agreement  with  the 
results  using  the  analytic  solutions  given  by  Roark.   This  is 
as  should  be  expected,  since  the  elastic  range  solution  using 
the  present  method  is  based  on  a  closed  form  analytic  solution. 

Results  in  the  plastic  range,  while  reasonable,  are  not  in 
as  close  agreement  with  other  published  results.   The  only  case 
where  data  for  a  direct  comparison  could  be  found  is  the 
configuration  discussed  in  Section  7.1.   From  Figure  7-3,  the 
deflections  at  loads  of  652  psi  and  852  psi  for  the  present 
solution  are  in  excellent  agreement  with  the  finite  element 
results  of  Armen  et  al. 

Results  for  the  annular  plate  of  Section  7.1  at  the  1052 
psi  load  are  not  as  encouraging.   Deflections  and  the  extent  of 
the  plastic  zone  are  greater  for  the  solution  obtained  by  Armen 
than  for  the  present  solution.   To  explain  this  difference,  the 
size  of  the  integration  and  load  increments  has  been  varied, 
but  it  shows  little  effect  on  the  results.   The  reason  behind 
this  difference  is  not  known.   However,  it  is  believed  that  it 
might  be  due  to  a  programming  error  either  in  this  report  or  in 
the  reference.   Also,  differences  could  arise  by  not  taking 
into  account  additional  integration  terms  arising  from  the 
singularity  at  r=r0,  as  discussed  in  Appendix  A. 
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Finally,  it  is  of  interest  to  compare  the  results  for  the 
simply  supported/guided  plate  of  Figure  7-9  to  the  results  for 
a  simply-supported  circular  plate  obtained  by  Moshaiov  and 
Vorus.   It  is  recognized  that  the  two  cases  are  not  equivalent. 
However,  one  would  expect  the  two  plates  to  behave  globally  in 
a  similar  fashion.   This  is  in  fact  the  case,  with  both 
deflections  and  plastic  zones  showing  reasonably  good  agreement 
in  view  of  the  differences  between  the  two  configurations. 
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CHAPTER  8 


SUMMARY  AND  CONCLUSIONS 


8 . 1  Summary 

This  thesis  has  developed  a  formulation  for  solving 
axisymmetr ically  loaded  annular  plate  bending  problems  using 
boundary  integrals.   This  particular  formulation  is  unique  in 
that  it  treats  the  annular  plate  as  a  one-dimensional  problem, 
using  "ring"  type  Green's  functions  to  determine  unknown 
boundary  conditions  and  arrive  at  a  plate  bending  solution. 
The  formulation  allows  a  numerical  incremental  load  method  to 
be  used  in  the  plastic  range  for  the  treatment  of  non-linear 
behavior.   Results  in  the  elastic  range  show  excellent 
agreement  with  results  obtained  using  a  conventional  analytic 
solution.   In  the  plastic  range,  results  are  in  reasonable 
agreement  with  those  obtained  using  a  finite  element  method. 

8 . 2  Conclus  ions 

This  thesis  treats  the  two-dimensional  problem  of  analysis 
of  axisymmetr ical ly  loaded  annular  plates  as  a  one-dimensional 
problem.   Using  a  method  of  solving  one-dimensional  problems 
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with  boundary  integrals  developed  by  Butterfield,  the  thesis 
demonstrates  that  a  closed  form  solution  for  the  annular  plate 
problem  can  be  obtained.   It  further  demonstrates  that  this 
approach  is  suitable  for  the  use  of  an  incremental  load  method 
to  obtain  a  solution  in  the  plastic  range. 

The  formulation  developed  in  this  thesis  is  advantageous 
over  a  conventional  analytic  solution  in  that  it  allows  a  wide 
range  of  problems  with  different  loading  and  boundary 
conditions  to  be  solved  using  a  single  algorithm.   Further,  the 
simplicity  of  the  approach  is  useful  from  an  educational 
standpoint  in  the  teaching  of  boundary  element  methods. 

8 .  3   Recommendations  for  Future  Work 

There  remains  room  for  much  additional  work  in  this  area. 
Some  suggestions  follow: 

1.  A  formulation  for  axisymmet rically  loaded  continuous 
circular  plates  remains  to  be  developed.   The  work  presented 
here  lays  a  foundation  for  the  continuous  plate  formulation. 
However,  extending  this  method  to  the  continuous  plate  case  may 
require  some  modifications. 

2.  In  developing  a  formulation  for  the  simple  beam 
problem  using  boundary  integrals,  Butterfield  uses  both  a 
direct  and  an  indirect  method,  and  shows  that  they  yield  the 
same  result.   Only  the  direct  method  is  presented  in  this 
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thesis.   The  indirect  method  is  a  more  intuitive  approach  than 
the  direct  method,  and  offers  advantages  in  furthering  the 
understanding  of  a  boundary  integral  solution.   It  would, 
therefore,  be  worthwhile  to  pursue  an  indirect  method  as  well 
as  a  direct  method  in  developing  a  boundary  integral 
formulation  for  circular  plate  geometries. 

3.  The  Green's  functions  selected  for  this  anlysis 
involve  many  terms  and  are  awkward  from  a  computational 
standpoint.   More  thought  needs  to  be  given  to  selecting  a 
"ring"  type  Green's  function  that  is  simpler  and  therefore 
offers  advantages  in  simplifying  calculations  and  reducing 
computational  time. 

4.  For  the  plastic  range,  the  lack  of  close  agreement 
between  the  results  presented  in  this  work  and  those  obtained 
by  Armen  et  al  using  a  finite  elemet  method  need  to  be  better 
understood.   Additional  comparative  data  should  be  found  or 
developed  to  increase  confidence  in  the  results  presented  in 
this  work. 


APPENDIX  A 


SELECTED  GREEN'S  FUNCTIONS 


This  Appendix  presents  the  selected  Green's  functions  and 
required  derivatives  that  are  used  to  solve  the  plate  bending 
problems  in  this  analysis.   The  functions  which  mathematically 
describe  the  deflection,  slope,  radial  and  tangential  moments, 
and  shear  are  obtained  from  Roark. 

A. 1   Description  of  Geometry  and  Sign  Conventions 

As  discussed  in  Chapter  4,  the  first  Green's  function 
selected  for  this  analysis  describes  the  response  of  annular 
plate  that  is  simply  supported  at  the  outer  radius  and 
unsupported  along  the  inner  radius  to  a  ring  load.   Figure  A-l 
depicts  this  plate  configuration  and  important  notation. 
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Figure  A-l.   Representation  of  Green's  function  used  in 
the  analysis. 


Roark  uses  a  sign  convention  with  deflection  wq  positive 
upward,  slope  #q  positive  when  the  deflection  wq  increases 
positively  as  r  increases,  moment  Mr  positive  when  creating 
compression  on  the  top  surface  of  the  plate,  and  the  shear 
force  Qqt  positive  when  acting  upward  on  the  inner  edge  of  an 
annular  section.   Subscripts  c  and  d  refer  to  the  radial 
pos  it  ion . 

This  sign  convention  differs  from  the  sign  convention  of 
Timoshenko,  which  has  been  adopted  for  this  analysis. 
Accordingly,  adjustments  must  be  made.   Specifically, 
Timoshenko  takes  deflection  to  be  positive  downward  and  shear 
to  be  positive  when  acting  downward  on  the  inner  edge  of  the 
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plate,  such  that  the  signs  for  these  items  as  given  by  Roark 
must  be  changed. 

For  slope,  Roark  uses  a  convention  opposite  to  that  of  Eq. 
(2-2-1)  and  the  convention  used  by  Timoshenko.   This 
difference,  combined  with  the  difference  in  defining  the  sign 
of  deflection,  effectively  cancel  each  other,  so  the  sign  of 
slope  from  Roark  does  not  change  for  this  analysis.   Finally, 
for  moments,  Roark  and  Timoshenko  use  the  same  sign  convention, 
so  no  sign  changes  for  moment  terms  are  necessary. 

The  formulas  in  the  sections  which  follow  include  a  number 
of  general  plate  functions  and  constants  that  are  not  depicted 
in  Figure  A-l.   Functions  L  and  G  and  their  derivatives  are 
included  in  Section  A. 7;  constants  F  and  C  are  included  in 
Sect  ion  A . 8 . 

A . 2   S  ingulari t  ies 

It  is  not  possible  to  evaluate  the  Green's  function  and 
associated  derivatives  at  the  exact  radial  location  where  the 
ring  load  is  applied,  due  to  singularities.   To  avoid  this 
problem  on  the  boundaries,  the  ring  load  is  positioned  a  small 
distance  inside  the  outer  plate  radius  and  outside  the  inner 
plate  radius  (0.00001  inches  for  the  results  presented  in 
Chapter  7).   During  final  review  of  this  thesis,  it  was  noted 
that  because  of  the  singularity  at  r=r0,  the  domain  integrand 
involving  plastic  moments  of  Eqs.  (5-4-2)  through  (5-4-6) 
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includes  singularities.   Therefore,  the  limiting  values  should 
be  explored  and  a  correction  applied  in  way  of  the  singular- 
ities.  Insufficient  time  was  available  to  repeat  the  analysis 
with  the  applied  correction  to  see  its  effect  on  the  results. 
Note  that  this  correction  affects  only  the  solution  in  the 
plastic  range. 

A . 3   Deflection.  Slope,  Moment  and  Shear 

The  following  expressions  describe  deflection  wq ,  slope 
$G,  radial  moment  Mgr,  and  shear  Qq    corresponding  to  the  first 
Green's  function: 
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A . 5   Second  Derivative  of  Deflection,  Slope,  Moment  and  Shear 

The  following  expressions  describe  the  second  derivatives 
with  respect  to  r0  of  wq,  $q,  Mgr  and  QQr 
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A. 6   Third  Derivative  of  Deflection,  Slope,  Moment  and  Shear 

The  following  expressions  describe  the  third  derivatives 
with  respect  to  r0  of  wq ,  #q,  Mqv    and  Qqc 


,3 

d    w. 


dr 


,3 

d    w 


dr 


Gd 


d3* 


dr 


Gd 


r    F 


r    d    G 
o 


d3* 


dr 


d3* 


dr 


Gd 


F         + 

4 


r2    d3G 


PG 


D       dr 


d3M 


dr 


Gr 


d3*rH    D 
—    Gd    -       F 

dr  r 

o 


PG    r 


d3G 


dr 


d3Q 


dr 


Gr 


88 


where 


,3 

d  w 


Gd 


dr 


D 


d3L 


1   A        3 

dr 
o 


d3L 


dr 


d3* 


dr 


Gd 


2      3 

p.c    d  L 


D  C_   dr 


A. 7   Plate  Functions  L3,  L9,  G3,  G6  and  G9  and  Their 
Der ivat  i ves 
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A.  8   Plate  Constants  F,.  F4.  F7.  d,  and  C7 
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A. 9   Derivatives  of  the  Green's  Function  With  Respect  to  r 

The  equations  describing  non-linear  plate  behavior 
included  in  Section  5-4  are  in  terms  of  the  first  and  second 
derivatives  of  the  Green's  function  wq  with  respect  to  r. 
These  derivatives  can  be  expressed  with  the  aid  of  Eqs.  (2-1-1) 
and  (2-1-7)  in  terms  of  moments  and  slopes  as  follows: 
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Expressions  for  slopes  and  moments  for  the  Green's 
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function  ( #q  and  MQr)  are  available  and  have  been  presented  in 
Section  A. 3.   It  is  further  noted  that  Eqs.  (5-4-2)  through  (5- 
4-6)  require  derivatives  of  the  above  two  expressions  with 
respect  to  r0.   This  is  accomplished  via  the  expressions  for 
derivatives  of  $q  and  Mqv    with  respect  to  r0  included  in 
Sections  A. 5  through  A. 7  above. 
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APPENDIX  B 
COMPUTER  PROGRAMS 


This  Appendix  contains  the  computer  programs  developed  to 
perform  the  analysis  work.   A  flow  chart  showing  the  overall 
program  structure  is  presented  in  Chapter  6. 

The  Appendix  is  organized  with  the  programs  INPUT  and  LOAD 
presented  first,  followed  by  the  program  MAIN.   The  supporting 
subroutines  for  the  program  MAIN  follow  that  program. 
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c 
c* 

C*** *t*tttt**t*tt*t****ttttttt*tt*t*t*tttt*tt***t****t*t* ************* 

c* 

C*  PROGRAM  INPUT 

C* 

C*  THIS  PROGRAM  ALLOWS  THE  USER  TO  DEFINE  THE  PROBLEM  OF  INTEREST 

C*  AND  TO  SET  CERTAIN  ADJUSTABLE  PARAMETERS  SUCH  AS  NUMBER  OF  DEPTH 

C*  INTEGRATION  INCREMENTS.    LOADING  CONDITIONS  ARE  INPUT  USING  THE 

C*  PROGRAM  LOAD.    TO  RUN,  THIS  PROGRAM  MUST  ACCESS  THE  FILES  GEOM.DAT, 

C»  BC.DAT.BCPOS.DAT  AND  ADJ. DAT 

C* 

Cttt****tt*t**t*tt**t*t**t**ttt*tttt*ttttttt*********t**t******t**t*** 

c* 
c 

PROGRAM  INPUT 
C 

IMPLICIT  REAL*8  (A-H.O-Z) 

CHARACTERS  BFLAG1  (  4  )  ,  BFLAG2  (4  )  ,  RESP 

CHARACTER*19  LBL ( 8) , LGEOM( 9 ) 

CHARACTER*25  LADJ(6) 
C 

REAL*8  ADJ(7) ,BC(4) ,GEOM(9) ,NU 
C 

INTEGER  POS(8) ,BCPOS(8) ,M,N 


C 
C* 

C*********************************#******************************** 

c* 

C*      INPUT  MATERIAL  CONSTANTS  AND  PLATE  GEOMETRY.    THE  FLEXURAL 
C*   RIGIDITY  IS  CALCULATED  BASED  ON  THE  INPUT  PARAMETERS.   N  IS  THE 
C*   TOTAL  NUMBER  OF  PARAMETERS 
C* 
C*******:M********************:M*********************************** 

c* 
c 

N=9 
C 

LGEOM(l)='OUTER  RADIUS  A  =  ' 

LGEOM(2)=' INNER  RADIUS  B  =  ' 

LGEOM(3)='POISSON  RATIO  NU  =  * 

LGE0M(4)='PLATE  CONSTANT  D  =  ' 

LGEOM(5)='THICKNESS  T  =  ' 

LGEOM(6)=* YIELD  STRESS  =  ' 

LGEOM(7)=' YIELD  STRAIN  =  ' 

LGEOM(8)='ULTIMATE  STRESS  =  ' 

LGE0M(9)='ULTIMATE  STRAIN  =  ' 
C 

WRITE(«, 1) 
1      F0RMAT(///,6X, 'WOULD  YOU  LIKE  TO  SEE  THE  CURRENT  PLATE  DIMENSIONS 
+     ,/,2X,*AND  MATERIAL  CONSTANTS  (Y/N)  ?  ' \) 

READ(*,30)  RESP 
C 

IF  (RESP  .EO.   'Y')  THEN 
WRITE(*,2) 
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2        FORMAT (//, 21, 'THE  CURBBNT  VALDBS  FOLLOW.   NOTE  THAT  A  ZERO*, 

+       /,2X,'VALDE  IS  SHOWN  FOB  D  WHEN  T  AND  E  ARE  GIVEN  AND  VICE', 
+       /,2X, 'VBHSA  ...',/) 

0PBN(6,FILE='GB0M.DAT' , STATUS =* OLD  * ) 
WRITE(*,30)  '  ' 
DO  8  1=1, N 

READ(6,70)  GBOM(I) 
WRITB(*,6)  LGBOM(I),GEOM(I) 
6  FORMAT(2X,A22,F20.10) 

8        CONTINUE 
CLOSE (6) 
ENDIF 


10 


12 


WRITE(*,9) 

F0RMAT(//,6X, 'DO  YOU  WISH  TO  INPUT  OR  CHANGE  GEOMETRY  OB', 
y         /,2X, 'MATERIAL  CONSTANTS  (Y/N)  ?  ',\) 
READ(*,30)  RBSP 

IF  (RESP  .BQ.  *Y')  THEN 

WRITE(*. 10) 

FORMAT(/,2X,  'INPUT  THE  VALUES  USING  DECIMAL  NOTATION  ',/) 

DO  12  1=1,3 

WRITE(*,20)  LGBOM(I) 

BBAD(*,70)  GEOM(I) 
CONTINUE 


13 


DO  13  1=5,9 

WRITB(*,20)  LGEOM(I) 

READ(*,70)  GEOM(I) 
CONTINUB 
GEOM(4)=GEOM(6)/GEOM(7)*GEOM(5)**3.0/12.0/(1.0-GEOM(3)**2.0) 


15 


0PEN(6,FILE='GE0M.DAT' , STATUS =' NBW* ) 
DO  15  1  =  1, N 

WRITE(6,70)  GEOM(I) 
CONTINUE 
CLOSE (6) 
ENDIF 


C 
C 

C **************************************************** ************** 

c 

C   INPUT  PLATB  BOUNDARY  CONDITIONS 
C 

C***»* ************************* ************************************ 

c 


LBL(1)='SHEAR  AT  A 
LBL(2)='MOMENT  AT  A 
LBL(3)='SLOPE  AT  A 
LBL(4)='DEFLBCTION  AT  A 
LBL(5)='SHEAR  AT  B 
LBL(6)='MOMENT  AT  B 
LBL(7)='SLOPE  AT  B 
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LBL(8)='DBFLBCTI0N    AT    B    =       ' 

OPEN ( 14, FILB=' BCLBL.DAT* , STATUS= ' NEW ) 
DO    155    1=1,8 

WHITB(14,156)    LBL(I) 

155  CONTINUE 

156  F0HMAT(A22) 
CLOSE(14) 

WRITE(*,160) 
160   FORMAT(//,6X, 'WOULD  YOU  LIKE  TO  SEE  THB  CURRENT  PLATE  BOUNDARY' 
+    ,/,2X, 'CONDITIONS  (Y/N)  ?  *\) 
READ(*,30)  RESP 

■ 

IF  (RESP  .EQ.  'Y')  THEN 

OPEN(ll,FILE=*BC.DAT' , STATUS= ' OLD ' ) 
OPBN ( 12, FILB=' BCPOS.DAT' , STATUS= ' OLD ' ) 
WRITE(*,30)  '  ' 
DO  165  1=1,4 

READ(11,70)  BC(I) 
165      CONTINUE 

DO  190  1=1,8 

READ(12,100)  BCPOS(I) 

IF  (BCPOS(I)  .GT.  0)  THEN 

WRITE(*,186)  LBL(I),BC(BCPOS(I)) 
186  FORMAT(2X,A22,F20.10) 

ENDIF 
190      CONTINUE 
CLOSE(ll) 
CLOSE(12) 
ENDIF 

WRITE(*, 17) 
17    FORMAT(//,6X, 'DO  YOU  WISH  TO  INPUT  OR  CHANGE  PLATE  BOUNDARY', 
+    /,2X, 'CONDITIONS  (Y/N)  ?  *,\) 
READ(*,30)  HBSP 

IF  (RESP  .EQ.  'Y')  THEN 
WRITB(«, 19) 
19    FORMAT(//,6X, 'INDICATE  WHETHER  THB  BOUNDARY  CONDITIONS  FOR  THE', 
€   /,2X,'PLATE  ARE  KNOWN  OR  UNKNOWN  BY  TYPING  IN  THB  LBTTER  "U*"  , 
8   /,2X,'FOR  UNKNOWN  BOUNDARY  CONDITIONS  AND  HITTING  RETURN  FOR*, 
e   /,2X, 'KNOWN  BOUNDARY  CONDITIONS',//) 

* 

WRITE(*,20)  'SHEAR  AT  A  ?   ' 
READ(*,30)  BFLAGl(l) 
WRITE(*,20)  'MOMENT  AT  A  ?   ' 
READ(*,30)  BFLAG1(2) 
WRITE(*,20)  'SLOPE  AT  A  ?   * 
READ(*,30)  BFLAGK3) 
WRITE(*,20)  'DEFLECTION  AT  A  ?   ' 
READ(*,30)  BFLAG1(4) 

WRITE(*,20)  'SHBAR  AT  B  ?   * 
READ(*,30)  BFLAG2(1) 
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WRITB(*,20)  'MOMENT  AT  B  ?   * 

READ(*,30)  BFLAG2(2) 

WRITE(*,20)  'SLOPE  AT  B  ?   ' 

BEAD(»,30)  BFLAG2(3) 

WRITE(*,20)  'DEFLECTION  AT  B  ?   * 

READ(*,30)  BFLAG2(4) 
20    FORMAT(2X, A22,\) 
30    FORMAT(Al) 

WRITE(*,40) 
40    F0RMAT(//,6X, 'TYPE  IN  THE  KNOWN  BOUNDARY  CONDITIONS  USING  ',/, 
«2X, 'DECIMAL  NOTATION  WHBN  PROMPTED ',// ) 


C 


c 


OPEN(8, F I LE=' BFLAG1.DAT' , STATUS =' NEW ' ) 
OPEN ( 9, FILB=' BFLAG2.DAT' , STATUS =' NEW ' ) 
OPEN(10, FILE='BC.DAT' , STATUS =' NEW  * ) 
C 

DO  50  1=1,4 

WRITE(8,30)  BFLAGl(I) 
IF  (BFLAGl(I)  .NE.  *U')  THEN 
L=L+1 

WHITE(*,20)  LBL(I) 
READ(*,70)  BC(L) 
WRITE(10,70)  BC(L) 
ENDIF 
50    CONTINUE 
C 

DO  60  1=1,4 

WRITE(9,30)  BFLAG2(I) 
IF  (BFLAG2(I)  . NE .  *U')  THEN 
L  =  L+1 

WRITB(*,20)  LBL(I+4) 
READ(*,70)  BC(L) 
WRITE(10,70)  BC(L) 
ENDIF 
60    CONTINUE 
C 

70    FORMAT(F20. 10) 
C 

CLOSE(8) 
CL0SE(9) 
CLOSE(IO) 
C 

OPEN(ll,FILE='POS.DAT' , STATUS =' NEW ' ) 
OPEN ( 1 2, F I LE=' BCPOS.DAT' , STATUS =' NEW ' ) 
C 

L  =  0 

DO  80  1=1,4 

IF  (BFLAGl(I)  .EQ.  'U')  THEN 
L=L+1 
POS(I)=L 
BCPOS(  11=0 
ELSE 
M  =  M+1 
BCPOS( I)=M 
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POS(I)=0 
BNDIF 
C 

WRITE(ll.lOO)  POS(I) 
WRITEC12. 100)  BCPOS(I) 
C 

80    CONTINUE 
C 

DO  90  1=5,8 

IF  (BFLAG2(I-4)  .EO.  *U')  THEN 
L=L+1 
POS(I)=L 
BCPOS(I)=0 
ELSE 
M  =  M+1 

BCPOS(I)=M 
POS(I)=0 
ENDIF 
C 

WRITE(ll.lOO)  POS(I) 
WRITE(12, 100)  BCPOS(I) 
C 

90   CONTINUE 
CLOSE(ll) 
CLOSEU2) 
100   FORMAT(Il) 
BNDIF 
C 

C* 
C******:M****************************************** **************** 

c* 

C»      INPUT  CALCULATION  ADJUSTMENT  ITEMS.   N  IS  THE  NUMBER  OF 

C*  INPUT  PARAMETERS  IN  THIS  SUBMODULE.   OBSERVE  THAT  SUBSEQUENT  TO 

C*  THE  DEVELOPMENT  OF  THIS  MODULE,  PARAMETER  NUMBER  3  WAS  NO  LONGER 

C*  ACCESSED  FROM  ANY  OTHER  PROGRAM,  AND  THEREFORE  SERVES  ONLY  AS  A 

C*  DUMMY  PARAMETBR.    THE  REMAINING  PARAMETERS  HAVE  THE  FOLLOWING 

C*  MEANINGS 

C* 

C*      1.   EPSILON  REFERS  TO  THE  DISTANCE  BETWEEN  THE  RING  LOAD 

C*  AND  THE  BOUNDARY.   NOTE  THAT  IN  SOME  INSTANCES,  LETTING 

C*  EPSILON  BQUAL  TO  ZERO  CAUSES  A  DIVIDE  BY  ZERO  TYPE  ERROR. 

C* 

C*       2.   THB  INCREMENTS  FOR  R  REFER  TO  THE  STATIONS  ACROSS 

C*  THE  PLATE.   THIS  APPLIES  TO  THE  NUMBER  OF  INTEGRATION 

C*  INCREMENTS  FOR  DISTRIBUTED  LOADS,  AS  WELL  AS  THE  NUMBER 

C*  OF  DATA  POINTS  FOR  THE  RESULTS. 

C* 

C*       3.    THIS  PARAMETER  IS  NO  LONGER  USED 

C* 

C*       4.    THIS  REFERS  TO  THE  EXTENT  BY  WHICH  THE  GREENS 

C*  FUNCTION  PLATE  GEOMETRY  OVERHANGS  THE  ACTUAL  PLATE. 

C*  FOR  EXAMPLE,  A  VALUE  OF  0.5  MEANS  THAT  THE  INNER  RADIUS  OF 

C*  THE  RING  LOADED  PLATE  DESCRIBED  BY  THE  GREEN'S  FUNCTION 

C*  OVERHANGS  THE  ACTUAL  PLATE  BY  0.5  UNITS. 

C* 
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C*       3.    THE  LOAD  INCREMENT  IS  THE  PERCENTAGE  OF  THE 

C»  ORIGINALLY  SPECIFIED  LOAD  MAGNITUDE  THAT  IS  APPLIED  DURING 

C»  EACH  LOAD  INCREMENT  AS  THE  PLATE  IS  LOADED  IN  THE  PLASTIC 

C*  RANGE 

C» 

C*      6.   THIS  IS  THE  NUMBER  OF  INTEGRATION  INCREMENTS  USED  TO 

C*  DETERMINE  THE  PLASTIC  MOMENT  ACROSS  THE  PLATE. 

C» 

C*      7.   EMS  IS  THE  ROOT  MEAN  SQUARE  OF  THE  PLASTIC  MOMENT 

C*  INCREMENTS  FOR  EVERY  POINT  IN  THE  PLATE  WHERE  THE  PLATE  IS 

C*  NO  LONGER  ELASTIC.   THIS  IS  COMPARED  WITH  THE  RMS  VALUE  FOR 

C*  THE  PREVIOUS  INCREMENT,  TO  SEE  WHETHER  THE  SOLUTION  HAS 

C*  SUFFICIENTLY  CONVERGED.   IX  CONVERGENCE  MEANS  THESE  NUMBERS 

C*  AGREE  TO  WITHIN  LESS  THAN  IX 

C* 

Ct*t***ttt*xttt**t**ttttttttttttt*****tttt*t*tt*ttt***t*tttttttttt*tt** 

C* 

c 

N  =  7 

LADJ(l)=*EPSILON  FOR  RO  =  ' 
LADJ(2)=' INCREMENTS  FOR  R  =  ' 
LADJ(3)='WIDTH  INTG  INCR  =  ' 
LADJ(4)='BOUNDS  FOR  A  AND  B  =  ' 
LADJ(5)='LOAD  INCRBMENT  =  ' 
LADJ(6)='DEPTH  INTG  INCR  =  ' 
LADJ(7)='X  RMS  CONVERGENCE  =  ' 
C 

WRITE(*, 102) 
102   FORMAT(//, 6X, 'WOULD  YOU  LIKE  TO  SEE  THE  CURRENT  ADJUSTABLE' 
+    ,/,2X,  'PARAMETER  SETTINGS?  '\) 
READ(*,30)  RESP 
IF  (RESP  .EQ.  *Y')  THEN 

OPEN( 16, FILE- 'ADJ. DAT* , STATUS ='  OLD  '  ) 
WRITE(*,30)  '  ' 
DO  108  1=1, N 

READ(16,70)  ADJ(I) 
WRITE(*,106)  LADJ(I) , ADJ(I) 
106  FORMAT(2X, A25, F20. 10) 

108      CONTINUE 
CL0SE(16) 
ENDIF 
C 

WRITB(*. 110) 
110   F0RMAT(//,6X, 'DO  YOU  WISH  TO  INPUT  OH  CHANGE  ANY  COMPUTATIONAL*, 
+    /,2X, 'ADJUSTMENT  ITEMS  (Y/N)  ?  ',\) 
READ(*,30)  RESP 
C 

IF  (RESP  .EQ.  ' Y' )  THEN 
C 

WRITE(«, 120) 
120       FORMAT(/,'    INPUT  THE  VALUES  USING  DECIMAL  NOTATION  ....',/) 
DO  125  1=1, N 

WRITE(*,126)  LADJ(I) 
READ(*,70)  ADJ(I) 
125      CONTINUE 


100 


126 


F0RMAT(2X, A25,\) 


0PEN(13,FILE=* ADJ.DAT* , STATUS =' NEW ) 
DO  130  1=1, N 

WRITE(13,70)  ADJ(I) 
130      CONTINUE 
CL0SE(13) 

ENDIF 


WRITE(*. 135) 
135   FORMAT(///,' 


YOU  ARE  LEAVING  THE  INPUT  MODULE',///) 


END 
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c 

C*t*******s**************************s****>**************************** 

c* 

C*  PROGRAM  LOAD 

C* 

C*      THIS  PROGRAM  ALLOWS  THE  USER  TO  INPUT  THE  DESIRED  LOADING 

C»  CONDITION  FOR  THE  PLATE.   UNIFORM  DISTRIBUTED  LOADS,  RAMPS,  AND 

C*  PARABOLICALLY  DISTRIBUTED  LOADS  ARE  PERMITTED,  AS  WELL  AS  RING 

C*  LOADS.   TO  RUN,  THIS  PROGRAM  MUST  ACCESS  THE  FILE  LOAD. DAT. 

C» 

C*      USING  THE  INPUT  DATA,  THE  PROGRAM  CALCULATES  THE  TOTAL  LOAD 

C*  PER  UNIT  LENGTH  THAT  ACTS  ON  A  RING  OF  THE  PLATE  SURFACE  OR  A 

C*  WIDTH  DETERMINED  BY  ADJUSTABLE  PARAMETER  #3  OF  THE  INPUT 

C*  PROGRAM 

C* 

C*****tttt*t**t*tt**tttttt*tttt*tt*t**t*****tt*t*tt*tt*tttttttt*t*ttt* 

C* 

c 

PROGRAM  LOAD 
C 

IMPLICIT  RBAL*8  (A-H.O-Z) 
CHARACTERS  RESP 
CHARACTER*30  LBL(4) 
REAL*8  LD(50) 
INTEGER  N 
C 

LBL(1)='MAGNITUDE  AT  OUTER  RADIUS  =  ' 
LBL(2)='DISTANCE  FROM  PLATE  CENTER  =  ! 
»   LBL(3)='MAGNITUDE  AT  INNER  RADIUS  =  ' 
LBL(4)='DISTANCE  FROM  PLATE  CENTER  =  ' 
C 

0PEN(6,FILE=" LOAD. DAT' , STATUS= ' OLD ' ) 
READ(6,90)  LD 
N=INT(6.0+2.0*LD(6) ) 
CLOSE(6) 
C 

WRITE(*. 10) 
10    FORMAT(//,'   DO  YOU  WISH  TO  SEE  THE  LOADING  CONDITIONS', 
+     /,  '   (Y/N)  ?  \\) 
READ(*,60)  RESP 
IF  (RESP  .EQ.  'Y')  THEN 
WHITE(*. 15) 
15        FORMAT(//,'   THE  EXISTING  LOADING  CONDITIONS  FOLLOW  ...',//) 
DO  40  1=1,4 

WHITE(*,30)  LBL(I),LD(I) 
30  FORMAT(2X,A30, F20. 10) 

40        CONTINUE 

WRITE(*,60)  '  ' 

IF  (LD(5)  .EC  1.0)  THEN 

WRITE(*,80'i  'LOADING  IS  LINEAR' 
ELSE 

WRITE(*,80(  'LOADING  IS  PARABOLIC 
ENDIF 

WRITE(*,60)  '  ' 
J  =  5 
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c 


DO  45  I=7,6+INT(LD(6) ) 
J  =  J+2 

WBITE(*,60)  '  ' 

WRITE(*,43)  'MAGNITUDE  OF  EING  LOAD', 1-6,  '  =  *,LD(J) 
WRITE(*,43)  'LOCATION  OF  BING  LOAD', 1-6,  '  =  ',LD(J+D 
43  FOHMAT(2X, A24, I3.A3.F20. 10) 

45  CONTINUE 
ENDIF 

WHITE(*,46) 

46  FORMAT(//,2X, 'DO  YOU  WISH  TO  CHANGE  THE  LOADING  CONDITIONS', 
+    /,  '   (Y/N)  ?  \\) 

READ(*,60)  RESP 

IF  (RESP  .EQ.  'Y')  THEN 

WRITE(*,50) 
50    FORMAT(//,2X, 'DO  YOU  WISH  TO  CHANGE  THE  DISTRIBUTED  LOAD', 
+   /,'     (Y/N)  ?  *\) 

READ(*,60)  RBSP 
60    FORMAT(Al) 

IF  (RESP  .EQ.  'Y')  THEN 

WRITE(*,70) 
70       FORMAT(//,2X, 'SPECIFY  THE  LOAD  PROFILE  USING  DECIMAL', 
+       /, '        NOTATION  ...  * ,//) 
DO  75  1=1,4 

WRITE(*,80)  LBL(I) 
READ(*,90)  LD(I) 
75       CONTINUE 
80       FORMAT(2X, A30,\) 
90        FORMAT(F20. 10) 

WRITE(*, 100) 
100      FORMAT(//,2X, ' IS  THE  LOAD  LINEAR  OR  PARABOLIC  (L/P)  ?  *,\) 
READ(*,60)  RBSP 
IF  (RESP  .EQ.  'L')  THEN 

LD(5)  =  1.0 
ELSE 

LD(5)  =  2.0 
ENDIF 
ENDIF 

WRITE(*. 110) 
110   FORMAT(//,2X, 'DO  YOU  WISH  TO  MODIFY  ANY  RING  LOADS  (Y/N)  ?  ',\) 
READ(*,60)  RESP 
IF  (RESP  .EQ.   'Y')  THEN 
WRITE(*, 120) 
120       FORMAT(//,2X,  'THE  NUMBER  OF  RING  LOADS  (INTEGER)  =  *\) 

READ(*,130,  N 
130      FORMAT(I3) 

LD(6)=REAL(N) 
K  =  5 

DO  200  1  =  1,  N 
K  =  K^2 


C 
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WRITB(*,60)  '  ' 

WRITE(*,210)  'MAGNITDDE  OF  RING  LOAD  * , I ,  '  =  ' 

READt*, 90   LD(K) 

WRITE(*,210)  'LOCATION  OF  HING  LOAD  ',1,  '=  ' 

READ(*,90)  LD(K+1) 
200      CONTINUE 
210      FORMAT(2X, A23, 13, A2,\) 
ENDIF 

OPEN(7,FILE='LOAD.DAT* , STATUS =  ' NEW  '  ) 
DO  219  1=1,50 
WRITE(7,90)  LD(I) 

219  CONTINUE 
ENDIF 

WRITE(*,220) 

220  FORMAT(//,6X, 'YOU  ARE  LEAVING  THE  LOADING  MODULE',///) 
END 
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c 

c* 

C***********«********xt**»***s*******s************t***************** 

c* 

C*       CALL  THE  SUBROUTINE  YLD,  WHICH  WILL  RETURN  A  LOAD  MATRIX 

C*   SUCH  THAT  YIELDING  HAS  JUST  STARTED  TO  OCCUR  IN  THE  PLATE 

C* 

C*******************i*********************************t********»*t** 

C* 

C 

CALL  YLD(KK, IMAX.LDMAT1) 
C 

C* 
C**s******************»******>************************************** 

c* 

C*       RBPEAT  THE  PROCESS,  AND  THEN  CALL  UPDATB  TO  REVISE  OUR  STATE 
C*   MATRIX  TO  REFLECT  THE  CONDITION  OF  THE  PLATE  AT  THE  ONSET  OF 
C*   YIELDING 
C* 

C*****************************************^************************ 

c* 

CALL  AMATR( BFLAG1 , BFLAG2 , BC , BCPOS , POS , LDMAT1 , LHS , AMAT) 

CALL  GAUSS(LHS, AMAT.X) 

CALL  BCM(POS,BC,X,BCMAT) 

CALL  SOLVER(LDMATl , BCMAT) 

CALL  UPDATE (STT1, IMAX , OUT, LDMAT1 , LL , LDMATO ) 

CALL  RESUL(M,STT1, COUNT, LDMATO) 
C 
C* 

C* 

C*        NOW,  INPUT  AN  INCREMENTAL  LOAD,  AND  COMPUTE  THE  RESULTING 

C*   MOMENTS, DEFLECTIONS  ETC,  USING  AS  AN  INPUT  THE  PLASTIC  MOMENTS 

C*   (IF  ANY)  AND  OTHER  DATA  FROM  THE  PREVIOUS  LOAD  CASE 

C* 

C*s****t************************************************************ 

c* 
c 

DO  110  1=1, NN 

LDMAT1(I)=ADJ(5)*LDMAT0(I) 
110   CONTINUE 
C 

116   CONTINUE 
C 

C       WRITE(*,130)  STT(IMAX,9) ,STT(IMAX, 10) 
C  130   FORMAT(5X, "MRP  IN  =  '.F20.10,'    MTP  IN  =  '.F20.10) 

CALL  AMA TR (BFL AG 1.BFLAG2.BC, BCPOS, POS, LDMAT1 , LHS, AMAT) 
CALL  GAUSS! LHS, AMAT, X) 
CALL  BCM(POS,BC, X, BCMAT) 
CALL  S0LVER(LDMAT1, BCMAT) 
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c 
c* 

C**«***x*s************************************************************* 
c* 

C«  GO  INTO  THE  SUBROUTINE  PLAST  TO  CALCULATE  THE  INCREMENTAL  PLASTIC 

C»  MOMENTS  AT  EACH  STATION  ACROSS  THE  PLATE  WHBRE  YIELDING  HAS 

C»  OCCURRED.   THE  PLASTIC  MOMENTS  THAT  ARE  RETURNED  BY  PLAST  ARE 

C*  SGUARBD  AND  ADDBD  TO  THE  TOTAL  IN  THE  VARIABLE  RMS.   RMS  IS 

C*  COMPARED  TO  THE  RMS  VALUE  FROM  THE  PREVIOUS  ITERATION  (SAVED  AS 

C*  RMS1) 

C* 

C********************************************************************** 

c* 
c 

RMS1=RMS 
RMS=0.0 

IZ=NINT(ADJ(6) ) 
DO  200  1=1, NN 

STT1(I,1)=STT(I,9) 

STT1(I,2)=STT(I,10) 

CALL  PLAST(I) 

IF  (SBT(I.IZ)  .GE.  GE0M(6))  THEN 

RMS=RMS+STT(I,9)**2.0+STT(I, 10)**2.0 
ENDIF 
200   CONTINUE 
C 

WRITE (*, 810)  RMS 
810   FORMAT (5X, 'RMS  =  *,F20.10) 
C 

C* 

C*     THIS  WAS  THROWN  IN  TO  AVOID  A  DIVIDE  BY  ZERO,  UNTIL  THE  PROCESS 
C*   HAS  BEBN  "PRIMED" 
C* 
C 

IF  (RMS  .EQ.  0.0)  THEN 

RMS=1.0 
BNDIF 
C 

IF  (ABS(SQRT(RMS)-SQRT(RMS1))/SQRT(RMS)*100.0  .GT.  ADJ(7))  THEN 
DO  815  1=1, NN 

STT(I,9)=(STT(I,9)+STTl(I,l))/2.0 
STT(I,10)=(STT(I,10)+STTl(I,2))/2.0 
815     CONTINUE 
GOTO  116 
BNDIF 
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c 

c* 

C*****************t******»*****************t»*t********************* 

c* 

C*       CALL  THE  SUBROUTINE  UPDATE,  WHICH  WILL  UPDATE  OUB  STATE 

C*  MATRIX  AND  OUR  LOADING  MATRIX  TO  RBFLECT  THE  CONDITION  OF  THE 

C*  PLATE  AT  THE  END  OF  THE  LOADING  INCREMENT.   THE  PARAMETER 

C*  COUNT  KEEPS  TRACK  OF  LOAD  INCREMENTS  FOR  PURPOSES  OF  DISPLAYING 

C*  RESULTS.   WHEN  COUNT=10,  RESULTS  ARE  DISPLAYED  AND  WRITTEN  TO 

C*  FILE.  (IE  BVERY  10  LOAD  INCREMENTS) 

C* 

Ctttttxtt**t*tttt*****ttttttttt* x**ttt**t*ttttt**t****tt*t**t**t*ttt 

c* 
c 

CALL  UPDATE (STT1 , IMAX , OUT , LDMAT1 , LL , LDMATO ) 

CALL  RESUL(M,STT1, COUNT, LDMATO) 

IF  (COUNT  .EQ.  10.0)  THEN 
COUNT=0.0 

ENDIF 

COUNT=COUNT+1.0 

M=M+1 

WRITB(*,117)  M 
117   FORMAT(5X, "LOAD  INCR  =  ',13) 
C 

C* 
Ct***t*t*t*t*tt*tttt*tttttt**t*ttt*t**tttt**t**tttt***tt***tt***tt** 

c* 

C*       CHECK  TO  SEE  IF  BXCEEDING  ULTIMATE  STRAIN.   IF  SO,  KICK 

C*  OUT  OF  THE  PROCEDURE.   IF  NOT,  WE  SHOULD  INCREMENT  ON  UP  AGAIN 

C*  BY  USING  AS  OUR  INPUT  THE  INCREMENTAL  LOAD  MATRIX. 
C* 

C» 

IF  (OUT  .NE.  1)  THEN 
GOTO  116 

ENDIF 
C 

CALL  RESUL(M,STT1, COUNT, LDMATO) 
C 

CLOSE(15) 
C 

END 
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c 
c* 

C**************s***************************************t************* 

c* 

C*  SUBROUTINE  INTEG 

C* 

C*  THIS  SUBROUTINE  IS  CALLED  BY  THE  PROGRAM  MAIN  TO  CREATE  THE 

C*  ARRAY  LDMAT.   THIS  ARRAY  CONTAINS  THB  VALUBS  OF  THE  TOTAL  LOAD 

C»  FOR  BACH  SEGMENT  ALONG  THE  PLATE  RADIUS.   IT  IS  CONSTRUCTED 

C*  BY  USING  THE  AVERAGE  OF  THE  TWO  BNDPOINTS  FOR  DISTRIBUTED 

C*  LOADS  AND  ADDING  TO  THIS  VALUE  THE  MAGNITUDE  OF  ANY  RING  LOADS 

C*  THAT  ARE  SITUATBD  WITHIN  THB  SBGMENT. 

C* 

c»*****************»************************»************************ 

C* 

c 

SUBROUTINE  INTEG ( LD , LDMAT) 
C 

IMPLICIT  REAL*8  (A-H.O-Z) 

REAL*8  LDMAT(50) ,LD(50) ,NU 

INTEGER  N 
$INCLUDB: 'COMMON. FOR* 
C 

A=GEOM(l) 

B=GEOM(2) 

NU=GEOM(3) 

D=GEOM(4) 
C 

K=NINT(ADJ(2)) 

DEL=(A-B)/ADJ(2) 

R  =  B 

N=NINT(6.0+2.0*LD(6)) 
C 

C* 

C»****»************************************************************** 
C* 

C*     CHECK  TO  SBE  IF  THERE  IS  A  DISTRIBUTED  LOAD.   IF  NOT,  PROCEED 
C*   TO  THE  SBCTION  FOR  RING  LOADS 
C* 
C****»************************************** ************************* 

c* 
c 

IF  (LD(1)  .NE.  0.0  .OR.  LD ( 3 )  . NE .  0.0)  THEN 
C 

R  =  B 

SLP=(LD(1)-LD(3) )/(LD(2)-LD(4)) 
DO  100  1=1, K 
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c 
c* 

C**«s**t*s****t************************************tt****t*********** 

c* 

C*     SBCTION  FOR  LINEARLY  DISTBIBUTED  LOADS 
C* 

C******t***tt*ttttt**********tttt**tt***t**tttttttttt*t***ttttt*t*t*t 

c* 
c 

IF  (LD(5)  .EQ.  1.0  )  THEN 
C 

B=H+DEL 
C 

IF  (R-LD(2)  .GT.  DEL/10.0)  THEN 
LDMAT(I)=0.0 
C 

ELSE  IF  (R  .GT.  LD(4))  THBN 

LDMAT(I)=(LP*(R-LD(4)-DEL/2.0)+LD(3) )*DEL 

ELSE 

LDMAT(I)=0.0 

ENDIF 


C 
C* 

c* 

C*      SECTION  FOR  PARABOLIC  LOADS   (Y=X**2) 
C* 

c* 
c 

ELSE  IF  (LD(5)  .EQ.  2.0  )  THEN 
C 

R=R+DEL 

P=(LD(2)-LD(4))**2.0/4.0/(LD(l)-LD(3)) 
IF  (LD(1)  .GT.  LD(3))  THEN 

Y=(R-DEL/2.0-LD(4) ) **2 . 0/4 . 0/P+LD ( 3 ) 
ELSE 

Y=-(LD(2)-R+DEL/2.0)**2.0/4.0/P+LD(4) 
ENDIF 
C 

IF  (R-LD(2)  .GT.  DEL/10.0)  THEN 
LDMAT(I)=0.0 
C 

ELSE  IF  (R  .GT.  LD(4))  THEN 
LDMAT( I)=Y*DEL 

ELSE 

LDMAT( I)=0. 0 

ENDIF 
ENDIF 
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c 

100      CONTINUB 
ENDIF 
C 

C* 

C**********:M*****************t ************************************** 
C* 

C*     SECTION  FOR  ADDING  IN  BING  LOADS 
C* 

C******************************************************************** 
C* 

c 

J  =  5 

IF  (N  .NE.  0)  THEN 

DO  300  I=7,6+NINT(LD(6) ) 
J  =  J+2 

L=INT((LD(J+1)-B)/DEL)+1 
LDMAT(L)=LDMAT(L)+LD( J) 
300      CONTINUE 
ENDIF 


RETURN 
END 


C 
C* 

C******************************************************************** 

c* 

C*    SUBROUTINE  AMATR 

C* 

C*      THIS  SUBROUTINE  IS  CALLED  BY  THE  PROGRAM  MAIN  TO  CREATE  THE 

C*    4X4  A  MATRIX,  WHICH  IS  SOLVBD  USING  GAUSS  ELIMINATION  TO 

C*    DBTERMINB  THE  UNKNOWN  BOUNDARY  CONDITIONS. 

C* 

C******************************************************************** 

C* 

c 

SUBROUTINE  AMATR ( BFLAG1 , BFLAG2 , BC , BCPOS , POS , LDMAT1 , LHS , AMAT) 
C 

IMPLICIT  REAL*8  (A-H.O-Z) 

CHARACTERS  BFLAGK4),  BFLAG2(4) 

REAL* 8  BC(4)  , LHS ( 4 ) ,  INC , NU , GRNMT1 ( 4 , 6 )  ,GRNMT2(4,6)  , 
+      AM AT (4, 4) , LDMAT(4) , LDTOT(4) ,LDMAT1(50) ,MPRT0T(4) ,MPTT0T(4) , 
+      STTAV(50,2) 

INTBGER  POS(8) , BCPOS(8) ,M,N 
C 

$INCLUDE: 'COMMON. FOR' 
C 

A=GEOM(l) 
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B=GK0M(2) 
NU=GE0M(3) 
D=GE0M(4) 
INC=ADJ(1) 
C 

DO  101  1=1,4 
LHS(I)=0.0 
IF  (I  .EQ.  1  .OR.  I  .EQ.  3)  THEN 

RO=A-INC 
ELSE 

RO=B+INC 
ENDIF 
C 

B  =  A 

CALL  GRNFNC(R,R0,NU,D,GRNMT1) 
R  =  B 

CALL  GRNFNC(R,RO,NU,D,GRNMT2) 
C 
C* 

c*****************»*********************************************»**«* 

C* 

C*      RBDUCE  THE  TERMS  ORIGINALLY  ON  THE  RHS  TO  A  4X4  MATRIX 

C*   COEFFICIENTS  (AMAT)  AND  A  4X1  MATRIX  OF  NON-COEFFICIENT 

C*   TERMS  (LHS) 

C* 

C****************»*»****** I****************************************** 

C* 

c 

IF  (I  .EQ.  1  .OR.  I  .EQ.  2)  THEN 
K  =  0 
L  =  0 
DO  90  J=l,4 

IF  (BFLAGl(J)  .NE.  'U')  THEN 
K  =  K+1 

LHS(I)=GRNMT1(J, 1 ) *BC ( K ) *A* ( - 1 . 0 ) ** ( J-l ) +LHS ( I ) 
ELSE 

L  =  L+1 

AMAT(I,L)=GRNMT1(J, 1 ) *A»(-1 . 0) **( J) 
ENDIF 
73  FORMAT(2X,F10.5) 

90  CONTINUE 


C 


DO  100  J=l,4 

IF  (BFLAG2(J)  . NB .  'U')  THEN 
K  =  K+1 

LHS(  I)-GRNMT2(  J,  1  )*BC(K)  *B*(-1  )»*J-tLHS(I) 
ELSE 

L  =  L+1 

AMAT(I,L)=GRNMT2( J, 1 ) *B* ( -1 ) ** ( J-l ) 
ENDIF 
100  CONTINUE 


ELSE 
K^O 
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L  =  0 

DO  110  J=l,4 

IF  (BFLAGl(J)  .NE.  'U')  THEN 
K  =  K+1 

LHS(I)=GRNMT1(J,2)*BC(K)*A*(-1.0)**(J-1)+LHS(I) 
ELSE 

L  =  L+1 

AMAT(I,L)=GRNMT1(J,2)*A*(-1.0)**(J) 
ENDIF 
110  CONTINUE 

C 

DO  115  J=l,4 

IF  (BFLAG2(J)  . NE .  *U')  THEN 
K  =  K+1 

LHS(I)=GHNMT2(J,2)*BC(K)*B*(-1)**(J)+LHS(I) 
ELSE 
L=L+1 

AMAT(I,L)=GRNMT2(J,2)*B*(-1)**(J-1) 
ENDIF 
115  CONTINUE 

ENDIF 
C 

C* 

C******************:M********* **************************************** 
C* 

C*     MOVE  THE  TERMS  ORIGINALLY  ON  THE  LHS  TO  THE  RHS ,  FORM  A  NEW  AMAT 
C*   MATRIX  AND  A  NEW  LHS  MATRIX  WHICH  CORRESPOND  TO  THE  A  AND  B  MATRICES 
C*   OF  THE  SYSTEM  AX=B 
C* 

C********************************************************************** 
C* 

IF  (BFLAG1(4)  . NE .  *U'  .AND.  I  .EQ.  1)  THEN 

LHS(1)=LHS(1)+R0*BC(BCP0S(4) ) 
ELSE  IF  (BFLAG1(4)  .EQ.  'U*  .AND.  I  .EQ.  1)  THEN 

AMAT(I,POS(4) )=AMAT(I,P0S(4) )-RO 
ENDIF 
C 

IF  (BFLAG2(4)  . NE .  *U'  .AND.  I  .EQ.  2)  THEN 

LHS(2)=LHS(2)+RO*BC(BCPOS(8) ) 
ELSE  IF  (BFLAG2(4)  .EQ.  'U'  .AND.  I  .EQ.  2)  THEN 

AMAT(I,POS(8) )=AMAT( I,POS(8) )-RO 
ENDIF 
C 

IF  (BFLAG1(3)  . NE .  'U'  .AND.  I  .EQ.  3)  THEN 

LHS(3)=LHS(3)-RO*BC(BCPOS(3) ) 
ENDIF 
IF  (BFLAGK4)  .  NE .   'U'  .AND.  I  .EQ.  3)  THEN 

LHS<  3)=LHS(3)+BC(BCPOSi4) ) 
ENDIF 
IF  (BFLAGK3)  .EQ.   'U'  .AND.  I  .EQ.  3)  THEN 

AMAT( I,POS(3) )=AMAT(I ,POS(3) )+RO 
ENDIF 
IF  (BFLAGK4)  .EQ.   'U'   .AND.  I  .EQ.  3)  THEN 

AMAT( I.POSf 4) )=AMAT( I , POS ( 4 ) J-1.0 
ENDIF 
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IF  (BFLAG2(3)  .ME.  'U*  .AND.  I  .EQ.  4)  THEN 

LHS(4)=LHS(4)-H0*BC(BCP0S(7) ) 
ENDIF 
IF  (BFLAG2(4)  . HE .  'U'  .AND.  I  .EQ.  4)  THEN 

LHS(4)=LHS(4)+BC(BCPOS(8) ) 
ENDIF 
IF  (BFLAG2(3)  .EQ.  'U'  .AND.  I  .BQ.  4)  THEN 

AMAT(I,P0S(7) )=AMAT(I,P0S(7))+R0 
ENDIF 
IF  (BFLAG2(4)  .BO.  *D*  .AND.  I  .EQ.  4)  THEN 

AMAT(I,POS(8) )=AMAT(I,POS(8) )-1.0 
ENDIF 
C 

101   CONTINUE 
C 

C* 
C************************S*************************************** 

c* 

C*      FINALLY,  ADJUST  THE  LHS  TO  ACCOUNT  FOR  THE  EXTBRNAL  LOAD 
C*   AND  PLASTICITY  BY  CALCULATING  AND  SUBTRACTING  APPROPRIATE 
C*   TBRMS. 
C* 

C******t**t***********************************************t****** 

c* 

c 

K=NINT(ADJ(2)) 
DBL=(A-B)/ADJ(2) 
C 

DO  500  1=1 ,K 

R=B-DEL/2.0+REAL(I)*DEL 
RO=A-ADJ(l) 

STTAV(I,1)=(STT(I,9)+STT(I+1,9) )/2.0 
STTAV(I,2)=(STT(I, 10)+STT(I+l,10))/2.0 
CALL  GRNFNC(R,R0,NU,D,GRNMT1) 

LDTOT( 1)=LDT0T(1)+LDMAT1(I)*GRNMT1(1, 1)*R 
MPRTOT(l)=MPRTOT(l)+STTAV(I,  1  )  *  ( -GRNMT1  (  3  ,  D/D  + 
+  NU/R*GRNMT1(2, 1))*DBL*R 

MPTTOT(l)=MPTTOT(l)+STTAV(I,2)*(-GRNMTl(2,  1)*DEL) 
LDT0T(3)=LDT0T(3)+LDMAT1(I)*GRNMT1(1,2)*R 
MPRTOT(3)=MPRTOT(3)+STTAV(I, 1 ) * ( -GRNMT1 ( 3 , 2 ) /D+ 
+  NU/R*GRNMT1(2,2) )*DEL*R 

MPTTOT(3)=MPTTOT(3)+STTAV(I,2)*(-GRNMTl(2,2)*DEL) 
RO=B+ADJ(l) 
CALL  GRNFNC(R,R0,NU,D,GRNMT1) 

LDT0T(2)=LDT0T(2)+LDMAT1(I)*GRNMT1( 1.  1)*R 
MPRTOT(2)=MPRTOT(2)+STTAV( I, 1 ) * ( -GRNMT1 ( 3 , 1 )/D+ 
+  NU/R*GRNMT1(2, 1))*DEL*R 

MPTTOT(2)=MPTTOT(2)+STTAV(  1 ,  2  )  *  ( -GRNMT1  (  2  ,  1)*DEL) 
LDTOTi4)=LDTOT(4  1+LDMATl( I)*GRNMT1( 1,2)*R 
MPHTOT(4)=MPRTOT(4)+STTAV(I, 1 ) » ( -GRNMT1 ( 3 , 2 ) /D+ 
+  NU/R*GRNMT1(2,2) )*DEL*R 

MPTTOT(4)=MPTTOT(4)+STTAV(I,2)*(-GRNMTl(2, 2)*DEL) 
500   CONTINUE 
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DO  105  1=1,4 

LHS(I)=LHS(I)-LDTOT(I)+MPHTOT(I)+MPTTOT(I) 
105   CONTINUE 
C 

DO  1000  1=1,4 
LDTOT(I)=0.0 
MPRTOT(I)=0.0 
MPTTOT(I)=0.0 
1000  CONTINUE 
C 

RETURN 
END 


C 
C» 

C***************************************************************** 

c* 

c*  SUBROUTINE  GRNFNC 

c* 

C*  THIS  SUBROUTINE  CALCULATES  THE  VALUES  OF  DEFLECTION,  SLOPE, 

C*  MOMENT  AND  SHEAR  AND  ASSOCIATED  DERIVATIVES  FOR  THE  GREEN'S 

C*  FUNCTION,  AND  RETURNS  THESE  VALUES  TO  THE  MAIN  PROGRAM.   TO 

C*  INDICATE  DERIVATIVES  OF  PARAMETERS,  THE  NUMBER  1,  2,  OR  3  IS 

C*  TACKED  ON  TO  THE  END  OF  THE  VARIABLE  NAME.   HENCE,  THE  FIRST 

C*  DERIVATIVE  OF  L9  WRT  RO  IS  LSI 

C* 

C* 

C*ttt**t**tt*tt*t**tt***t*t*****t****ttt****t*tt*t*t***tttt**t**tttt 

c* 
c 

SUBROUTINE  GRNFNC ( R , HO , NU , D , GRNMAT) 
C 

IMPLICIT  REAL*8  (A-H.O-Z) 

RBAL*8  NU, GRNMAT (4, 6) , L3 , L31 , L32 , L33 , L9 , L91 ,  L92 , L93 
tINCLUDE: 'COMMON. FOR' 
C 

11=1.0 

A=GE0M(1)+ADJ(4) 

B=GEOM(2)-ADJ(4) 

THB=0.0 

MB  =  0.0 

QB  =  0.  0 
C 

IF  (R  .GT.  RO)  THEN 
RRO=1.0 

ELSE 

RRO  =  0.0 

ENDIF 
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c 
c* 

C**tt********************X******************************************* 

c* 

C*      CALCULATB  PLATE  CONSTANTS  AND  PLATE  FUNCTIONS,  AND  THEIR 

C*    DEBIVATIVBS 

C* 

C******************************************************************** 
C* 

c 

L3=RO/4.0/A*(((RO/A)**2.0+1.0)*DLOG(A/RO)*(RO/A)**2.0-1.0) 

L31=1.0/4.0/A*(DLOG(A/RO)*(3.0*RO*»2.0/A**2.0+1.0)+2.0* 
+       (fiO**2.0/A»*2. 0-1.0)) 

L32=1.0/4.0/A*(6.0*RO/A«*2.0*DLOG(A/RO)+RO/A**2.0-1.0/RO) 

L33=1.0/4.0/A*(6.0/A**2.0*(DLOG(A/RO)-1.0)+1.0/A**2.0+1.0/RO 
+      **2.0) 
C 

L9=RO/A*((1.0+NU)/2.0*DLOG(A/RO)+(1.0-NU)/4.0*(1.0-(RO/A)**2.0)) 

L91=1.0/4.0/A*(2.0*(1.0+NU)*DLOG(A/RO)-(3.0*NU+1.0)-3.0*(1.0-NU) 
+       *HO**2.0/A**2.0) 

L92=-(1.0+NU)/2.0/A/RO-3.0*RO*(1.0-NU)/2.0/A**3.0 

L93=(1.0+NU)/2.0/A/RO**2.0-3.0*(1.0-NU)/2.0/A**3.0 
C 

C1=(1.0+NU)/2.0*B/A*DLOG(A/B)+(1.0-NU)/4.0*(A/B-B/A) 

C7=1.0/2.0*(1.0-NU**2.0)*(A/B-B/A) 
C 

G3=RO/4.0/R*( ( (RO/H)**2.0+1.0)*DLOG(R/RO)+(RO/H)**2.0-1.0)*RRO 

G31=1.0/4.0/R*(3. 0*RO**2.0/R**2.0*DLOG(R/RO)+2.0*RO**2.0/R«*2.0+ 
+       DLOG(R/RO)-2.0)*RRO 

G32=1.0/4.0/R*(6. 0*RO/R**2 . 0*DLOG ( R/RO) +R0/R**2 . 0-1. 0/RO)*RRO 

G33=1.0/4.0/R*(6.0/R**2.0*(DLOG(R/RO)-1.0)+1.0/R**2.0+1.0/RO**2.0) 
+       *RRO 
C 

G6=RO/4.0/R*( ( RO/ R ) **2. 0-1. 0+2. 0*DLOG( R/RO) )*RR0 

G61=1.0/4.0/R*(3.0*(RO**2.0/R**2.0-1.0)+2.0*DLOG(R/RO) )*RRO 

G62=(3.0*RO/2.0/R**3.0-1.0/2.0/RO/H)»RRO 

G63=(3.0/2.0/R**3.0+1.0/2.0/R/RO**2.0)*RRO 
C 

G9=HO/R*((1.0+NU)/2.0*DLOG(R/RO)+(1.0-NU)/4.0*(1.0-(RO/R)**2.0)) 
+       *RRO 

G91 = 1. 0/4. 0/B* ( 2. 0*(1.0+NU)*(DLOG( R/RO )-1.0)+( 1.0-NU) -3.0* 
+       (1.0-NU)*RO**2.0/R**2.0)*RRO 

G92=-l. 0/4. 0/R*( 2. 0*(1.0+NU)/RO+6.0*RO/R**2.0*( 1.0-NU) )*RRO 

G93=1.0/4.0/R*(2.0*(1.0+NU)/RO**2.0-6.0/R**2.0*( 1.0-NU) )*RR0 
C 

F1=(1.0+NU)/2.0*B/R*DLOG(R/B)+(1.0-NU)/4.0*(R/B-B/R) 

F2=1.0/4.0*(1.0-(B/R)**2.0*(1.0+2.0*DLOG(R/B) ) ) 

F3=B/4.0/R*( ( (B/R)**2.0+1.0)*DLOG(R/B)+(B/R)**2. 0-1.0) 

F4=1.0/2.0*( (1.0+NU)*B/R+( 1.0-NU) *R/B) 

F5=1.0/2.0*(1.0-(B/R)**2.0) 

F6=B/4.0/R*( (B/R)**2.0-1.0+2.0*DLOG(R/B) ) 

F7=1.0/2.0*(1.0-NU**2.0)*(R/B-B/R) 

F8=1.0/2.0*( 1.0+NU+(1.0-NU)*(B/R)**2.0) 

F9=B/R*f (1.0+NU)/2.0*DLOG(R/B)-(1.0-NU)/4.0*( 1.0-(B/R)**2.0)) 
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c 

c» 

C***t************************»***************s*********************** 

c* 

C*  CALCULATE  THE  DEFLECTION,  SLOPE,  MOMENT  AND  SHEAR,  AND 

C*  ASSOCIATED  DERIVATIVES,  AND  STORE  IN  THE  MATRIX  GRNMAT. 

C*  THIS  MATRIX  IS  A  4  BY  4  MATRIX,  WHERE  ROW  POSITIONS  ARE  AS 

C*  FOLLOWS 

C* 

C*  1    DEFLECTION 

C*  2    SLOPE 

C*  3    MOMENT 

C*  4    SHEAR 

C* 

C*  AND  COLUMN  POSITIONS  ARE  AS  FOLLOWS 

C* 

C*  1    GREBN'S  FUNCTION  (NO  DERIVATIVES) 

C*  2    1ST  DERIVATIVE  OF  GREEN'S  FUNCTION  WRT  RO 

C*  3    2ND  DERIVATIVE  OF  GREEN'S  FUNCTION  WRT  RO 

C»  4    3RD  DERIVATIVE  OF  GREEN'S  FUNCTION  WRT  RO 

C* 

C*  FOR  EXAMPLE,  THE  VALUE  RETURNED  BY  GRNMAT(2,3)  REPRESENTS 

C*  THE  SECOND  DERIVATIVE  OF  THE  SLOPE  FOR  THE  GREEN'S  FUNCTION 

C*  WRT  RO 

C* 

c******************************************************************** 

c* 

c 

YB--W*A**3.0/D*(C1*L9/C7-L3) 

THB=W*A**2.0/D/C7*L9 

GRNMAT (1,1)=-(YB+THB*R*F1-W*R**3.0/D*G3) 

GRNMAT(2, 1 ) =THB*F4-W*R**2 . 0/D*G6 

GRNMAT ( 3, 1)=THB*D/R*F7-W*R*G9 

GRNMAT(4, 1 ) =W*RO/R*RRO 


YB1=-W*A**3.0/D*(C1*L91/C7-L31) 

THB1=W*A**2.0/D/C7*L91 

GRNMAT (1, 2 ) =- ( YB1+THB 1*R*F1-W*R**3 . 0/D*G3 1 ) 

GRNMAT(2,2)=THB1*F4-W*R**2.0/D*G61 

GRNMAT(3,2)=THB1*D/R*F7-W*R*G91 

GRNMAT(4,2)=W/R*RRO 


YB2=-W*A**3.0/D*(C1/C7*L92-L32) 

THB2=W*A**2.0/D/C7*L92 

GRNMAT ( 1,3)=-(YB2+R*F1*THB2-W*R**3. 0/D*G32) 

GRNMAT(2,3)=F4»THB2-W*R»*2.0/D*G62 

GRNMAT(3,3)=D*F7/R*THB2-W*R*G92 

GRNMAT(4,3)=0 


YB3=-W»A»*3. 0/D*vCl/C7*L93-L33) 

THB3-W*A**2.0/D/C7*L93 

GRNMAT(1,4)=-(YB3^R*F1*THB3-W*R**3.0/D*G33; 


118 


GRNMAT(2,4)=F4*THB3-W*R**2.0/D*G63 

GRNMAT(3,4)=D*F7/R*THB3-W*R*G93 

GRNMAT(4,4)=0 

HETDHN 
END 


C 

C* 

Ct************:M********** ******************************************* 

C* 

C*    SUBROUTINE  GAUSS 

C* 

C*      THIS  SUBROUTINE  USES  GAUSS  ELIMINATION  TO  DETERMINE  THE 

C*    UNKNOWN  BOUNDARY  CONDITIONS.   THE  MATRICES  A  AND  B  OF  THE 

C*    BQUATION  AX  =  B  ARE  PROVIDED  BY  THE  SUBROUTINE  AMATR.   THIS 

C*    SUBROUTINE  THEN  RETURNS  THE  MATRIX  X  OF  UNKNOWN  BOUNDARY 

C*    CONDITIONS.   THIS  SUBROUTINE  WAS  TAKEN  FROM  THE  BOOK 

C*    "NUMERICAL  ANALYSIS"  BY  L.W.  JOHNSON  AND  R.D.  RIESS. 

C* 

C******************************************************************** 

c» 
c 

SUBROUTINE  GAUSS(B,A,X) 
C 

IMPLICIT  REAL*8  (A-H.O-Z) 

CHARACTER*22  BCLBL(8) 

REAL*8  A(4,4) ,B(4) ,X(4) ,BCMAT(8) ,BC(4) 

DIMENSION  AUG(50,51) 

INTEGER  POS(8) ,BCPOS(8) ,M,N 
C 

NM1  =  3 

NP  =  5 

N  =  4 
C 

C   SET  UP  THE  AUGMENTED  MATRIX  FOR  AX=B 
C 

DO  2  1  =  1, N 
DO  1  J=1,N 
AUG(I, J)=A(I, J) 

1  CONTINUE 
AUG( I,NP)=B(I) 

2  CONTINUE 
C 

C   THE  OUTER  LOOP  USES  ELEMENTARY  ROW  OPERATIONS  TO  TRANSFORM 
C   THE  AUGMENTED  MATRIX  TO  ECHELON  FORM 
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DO  8  1=1, NM1 
C 

C   SEARCH  FOB  THE  LABGEST  ENTBY  IN  COLUMN  I,  BOWS  I  THHODGH  N 
C   IPIVOT  IS  THB  ROW  INDEX  OF  THE  LARGEST  ENTBY 
C 

PIVOT=0.0 

DO  3  J=I,N 

TKMP=ABS(AUG(J,I) ) 

IF(PIVOT.GE.TEMP)  GO  TO  3 

PIVOT=TEMP 

IPIVOT=J 

3  CONTINUE 
IF(PIVOT.BQ.O.O)  GO  TO  13 
IFdPIVOT.BQ.  I)  GO  TO  5 

C 

C   INTERCHANGE  ROW  I  AND  ROW  IPIVOT 

C 

DO  4  K=I,NP 

TEMP=AUG(I,K) 

AUG(I,K)=AUG(IPIVOT,K) 

AUG(IPIVOT,K)=TEMP 

4  CONTINUE 
C 

C   ZBRO  ENTRIES  ( 1+1 ),( i+2 ),...,( N , 1 )  IN  THE  AUGMENTED  MATRIX 
C 

5  IP1=I+1 

DO  7  K=IP1,N 
Q=-AUG(K, I)/AUG(I, I) 
AUG(K, I)=0.0 

DO  6  J=IP1,NP 

AUG(K, J)=Q*AUG(I, J)+AUG(K, J) 

6  CONTINUE 

7  CONTINUE 

8  CONTINUE 

IF(AUG(N,N) .EQ.0.0)  GO  TO  13 
C 

C   BACKSOLVE  TO  OBTAIN  A  SOLUTION  TO  AX=B 
C 

X(N)=AUG(N,NP)/AUG(N,N) 
DO  10  K=1,NM1 
Q=0.0 

DO  9  J=1,K 
Q=Q+AUG(N-K,NP-J)*X(NP-J) 

9  CONTINUE 
X(N-K)=(AUG(N-K,NP)-Q)/AUG(N-K,N-K) 

10  CONTINUE 
C 

C   CALCULATE  THE  NORM  OF  THE  RESIDUAL  VECTOR,  B-AX. 

C   SET  IERROR-1  AND  RETURN 

C 

RSO=0.0 

DO  12  1  =  1,  N 
Q=0.0 

DO  11  J=1,N 


120 


Q=Q+A(I,J)*X(J) 

11  CONTINUE 
HESI=B(I)-Q 
RMAG=ABS(RBSI) 
RSQ=RSQ+RMAG**2 

12  CONTINUE 
RNOBM=SQBT(BSQ) 
IERROR=l 

C 

C   ABNORMAL  RBTURN  REDUCTION  TO  ECHELON  FORM  PRODUCES  A  ZERO 

C   ENTRY  ON  THE  DIAGONAL.   THE  MATBIX  A  MAY  BE  SINGULAR. 
C 

13  IEBROR=2 
C 

RETUBN 
END 


C 
C* 

C******** ************************************************  ************ 

c* 

C*  SUBROUTINE  BCM 

C* 

C*      THIS  SUBROUTINE  TAKES  THE  KNOWN  BOUNDARY  CONDITIONS  STORED 

C*  BY  THE  PROGRAM  INPUT  IN  THE  FILE  BCM. DAT  AS  WELL  AS  THE 

C*  "UNKNOWN"  BOUNDARY  CONDITIONS  RETURNED  BY  THE  SUBROUTINE  GAUSS 

C*  AND  CREATES  THB  8X1  MATRIX  OF  BOUNDARY  CONDITIONS  BCMAT. 

C* 

C***** ******************************************************  ********* 

c* 


SUBROUTINE  BCM( POS , BC , X , BCMAT ) 

IMPLICIT  REAL*8  (A-H.O-Z)C 
REAL»8  BC(4) ,X(4) ,BCMAT(8) 
INTEGER  POS(8) 

K  =  0 
L  =  0 
DO  120  1=1,8 

IF  (POS(I)  .GT.  0)  THEN 
K  =  K+1 

BCMAT(I)=X(K) 
ELSE 

L=L+1 

BCMAT(I)=BC(L) 
ENDIF 
120   CONTINUE 
RETURN 
END 
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c 
c* 

C**************************************t***************************** 

c* 

C*    SUBROUTINE  SOLVER 

C* 

C*      THIS  SUBROUTINE  SOLVES  TAKES  THE  COMPLETE  SET  OF  BOUNDARY 

C*    CONDITIONS  SUPPLIED  BY  THE  SUBROUTINE  BCM  AND  SOLVES  THE  PLATE 

C*    BENDING  PROBLEM  ACROSS  THE  PLATE  WIDTH.   IMPORTANT  PARAMETERS 

C*    SUCH  AS  DEFLBCTION,  SLOPE  ETC  ARE  STORED  IN  THE  ARRAY  STT. 

C* 

C******************************************************************** 

C» 

C 

c 
c 


SUBROUTINE    SOLVER ( XMT , BCMAT) 
IMPLICIT    REAL*8    (A-H.O-Z) 


REAL*8    GRNMTH4.6)  ,GRNMT2(4,6)  ,BCMAT(8)  ,  DT  (  50  ,  4  )  ,  SM(  50  ,  4  )  ,  NU  , 
+       XMT(50) ,STTAV(50,2) 
INTEGER    M,N 
(INCLUDE: 'COMMON. FOR" 
C 

C* 
C«************»***************************************************** 

c* 

C*   REZERO  THE  MATRICES  DT  AND  ST 

C» 

C***t**ttt*t****t**t**tt***ttttt*tttttttttt**t*t**t**tttt*t**tt***tt 

C* 


DO  10  1=1,50 

DO  5  J=l,4 

SM(I,J)=0.0 

DT(I, J)=0.0 

5 

CONTINUE 

10 

CONTINUE 

C 

A=GBOM(l) 

B=GEOM(2) 

NU=GEOM(3) 

D=GE0M(4) 

C 

C 

C 

RO  =  B 

DEL= (A-B)/ADJ(2) 
N=INT(ADJ(2) ) 
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c 
c« 

C**t*********«************»*****t*s****«****«*t**tt******tt***x****t 

c* 

C»      THIS  LOOP  CALCULATES  THE  LOAD  TERM  AND  THE  PLASTIC  MOMENT 

C*  TERM  FOR  EACH  BQUATION.   I  REFERS  TO  THE  RADIAL  STATION  AND  J 

C»  REFERS  TO  THE  EQUATION,  WHERE  THE  FIRST  EQUATION  IS  FOR 

C*  DEFLECTION,  THE  SECOND  FOR  DW/DR,  THE  THIRD  FOR  D~2W/DR~2  AND 

C*  THE  FOURTH  FOR  D~3W/DH~3.   THE  DERIVATIVES  ARE  THEN  USED  TO 

C*  CALCULATE  SLOPBS,  MOMENTS  AND  SHEARS 

C* 

CS**********t************************************?***********»****** 

C* 

c 

DO  155  1=1, N+l 
C 

IF  (I  .EQ.  1)  THEN 

RO=B+ADJ(l) 
ENDIF 
IF  (I  .EQ.  (N+l))  THEN 

RO=A-ADJ(l) 
ENDIF 
DO  150  M=1,NINT(ADJ(2)) 

R=B-DEL/2.0+REAL(M)»DEL 
CALL  GRNFNC(R,R0,NU,D,GRNMT1) 
STTAV(M, 1)=(STT(M,9)+STT(M+1,9) )/2.0 
STTAV(M,2)=(STT(M, 10)+STT(M+l,10))/2.0 
DO  140  J=l,4 

DT( I, J)=DT(I, J)+(XMT(M)*GRNMT1(1, J  )  *R-DEL*STTAV (M , 1)* 
+  (-GRNMTK3,  J)/D  +  NL'/R*GRNMT1(2,  J)  )*R 

+  +DBL*STTAV(M,2)*GRNMT1(2,J) ) 

140         CONTINUE 
150      CONTINUE 

IF  (I  .EQ.  1)  THEN 

RO=B+DEL 
ELSE 

RO=RO+DEL 
ENDIF 
155   CONTINUE 
C 
C 
C* 

c* 

C*     AT  EACH  INCREMENT  OF  WIDTH  DEL  ALONG  THE  PLATE  RADIUS,  CALCULATE 

C»   THE  DEFLECTION,  SLOPE,  MOMENT  AND  SHBAR,  WHICH  ARE  STORED  IN 

C*   STT  (BACH  COLUMN  REPRESENTS  RO ,  Y,  TH,  M,  Q,  Y*  '  ,  Y'"'  CORRES- 

C*   PONDING  TO  EACH  STATION  I  ACROSS  THE  RADIUS! 

C* 

C********»*t***************»**»*************»************************* 

C* 

C 

DO  200  1=1, N+l 
C 

IF  fl  .EQ.  1)  THEN 
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B0=B+ADJ(1) 
BNDIF 
IF  (I  .BQ.  (N+l))  THEN 

B0=A-ADJ(1) 
BNDIF 
C 

B  =  A 

CALL  GBNFNC(B,B0,NU,D,GRNMT1) 

B=B 

CALL  GRNFNC(RtB0,NU,D,GRNMT2) 

DO  180  J=l,4 

DO  170  K=l,4 

SM(I,J)=A*GRNMT1(K,J)*BCMAT(K)*(-1.0)**K- 
+  B*GRNMT2(K, J ) *BCMAT ( K+4 ) * ( - 1 . 0 ) **K+SM( I , J) 

170         CONTINUB 

SM(I, J)=SM(I,J)+DT(I,J) 
180      CONTINUE 
C 
C* 

c* 

C*  STT(I.l)  -  BADIUS  OF  INTBRKST 

C*  STT(I,2)  -  DEFLBCTION 

C*  STT(I,3)  -  DW/DB 

C*  STT(I,4)  -  BADIAL  MOMENT 

C*  STT(I,5)  -  SHEAR 

C*  STT(I,6)  -  TANGENTIAL  MOMENT 

C*  STT(I,7)  -  DW"2/DR"2 

C*  STT(I,8)  -  DW*3/DRA3 

C*  STT(I,9)  -  BADIAL  PLASTIC  MOMENT 

C*  STT(I.IO)-  TANGENTIAL  PLASTIC  MOMENT 

C* 

C*tt*tt**t**t*t*t***t**t*t**********t*******t**t*t**t*****tt*t*****t* 

c* 
c 

STT(I,l)=BO 

STT(I,2)=SM(I,l)/RO 

STT(I,3)=1.0/RO*(SM(I,2)-STT(I,2)) 

STT(I,7)=1.0/BO*(SM(I,3)-2.0*STT(I,3)) 

STT(I,8)=1.0/BO*(SM(I,4)-3.0*STT(I,7)) 

STT(I,4)=-D*(STT(I,7)+NU/BO*STT(I,3)) 

STT(I,5)=D*(STT(I,8)+1.0/BO*STT(I, 7)-1.0/BO**2.0* 
+       STT(I,3)) 

STT(I,6)=-D*(1.0/BO*STT(I,3)+NU*STT(I,7)) 
C 

C* 

C*       THIS  IS  A  STEP  TO  PBEVENT  THE  INTERIOR  POINTS  FROM 
C*   HAVING  THE  EPSILON  THROWN  IN. 
C* 
C 

IF  (I  .EQ.  1)  THEN 
RO=B+DEL 

BLSE 

RO=RO+DEL 

ENDIF 
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c 

200   CONTINUE 
C 

RETURN 
END 


C 
C* 

C*x******************************s*********************************** 
c* 

C*  SUBROUTINE  YLD 

C* 

C*  THIS  SUBROUTINE  TAKES  THE  PLATE  BENDING  SOLUTION  FOR  THE  INPUT 

C*  LOAD  INCREMENT  AND  CALCULATES  THE  LOAD  AT  WHICH  YIELDING  OF  THE 

C*  OUTER  FIBER  WILL  OCCUR.   THE  LOAD  MATRIX  IS  THEN  INCREASED  TO 

C*  THIS  VALUE.   THIS  IS  ESSENTIALLY  A  TIME  SAVER,  TO  AVOID  STEPPING 

C*  UP  THE  LOAD  IN  THE  ELASTIC  RANGE. 

C* 

c* 
c 

c 


SUBROUTINE  YLD ( KK , IMAX , LDMAT1 ) 


IMPLICIT  REAL*8  (A-H.O-Z) 
HEAL*8  LDMATK50)  ,  STRMAX  ,  STTMAX 
INTEGER  KK, IMAX, IRMAX, ITMAX 
C 

SINCLUDE: 'COMMON. FOR' 
C 

C* 

C*  FIND  THE  LOCATION  ACROSS  THE  PLATE  WHERE  THE  BENDING  MOMENTS 

C*   ARE  THE  HIGHEST 
C* 
C 

STRMAX=0.0 
STTMAX=0.0 
DO  122  I=1,NINT(ADJ(2) )+l 

IF  (ABS(STT( 1,4) )  .GT.  STRMAX)  THEN 
STRMAX=ABS(STT(I,4) ) 
IRMAX=I 
ENDIF 

IF  (ABSCSTT-;  1,6)  )  .GT.  STTMAX'i  THEN 
STTMAX  =  ABS(STT(  1,6)) 
ITMAX=I 
ENDIF 

IF  (STRMAX  .GT.  STTMAX)  THEN 
KK=1 

IMAX=IRMAX 
ELSE 
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11=2 

IMAX=ITMAX 
KNDIF 
122      CONTINUE 
C 

C* 

C*     CALCULATE  THE  EQUIVALENT  STRESS  AT  THE  LOCATION  WHERE  STRESSES 
C»   ARE  HIGHEST.   COMPARE  IT  TO  THE  YIELD  STRESS.   THEN  RATIO  UP  THE 
C»   LOAD  MATRIX  ACCORDINGLY 
C* 
C 

ERMAX=STT(IMAX,7)*GEOM(5)/2.0 
ETMAX=STT(IMAX,3)*GEOM(5)/2.0/STT(IMAX, 1) 

SRMAX=GEOM(6)/GEOM(7)/(1.0-GEOM(3)**2.0)*(ERMAX+GEOM(3)*ETMAX) 
STMAX=GEOM(6)/GEOM(7)/(1.0-GEOM(3)**2.0)*(ETMAX+GEOM(3)*ERMAX ) 
SEMAX=SQRT(SRMAX**2.0-SRMAX*STMAX+STMAX**2.0) 
R=GBOM(6)/SEMAX 
DO  126  J=1,NINT(ADJ(2)) 
LDMAT1(J)=R*LDMAT1(J) 
126      CONTINUE 

WRITB(*,130)  LDMATl(l) 
130      FOHMAT(2X, 'LDMAT1  =  '.F20.10) 
C 

RETURN 
END 


C 
C* 

c* 

C*  SUBROUTINE  PLAST 

C* 

C*  THIS  SUBROUTINE  CALCULATES  PLASTIC  MOMENTS  USING  A  TRAPAZOIDAL 

C*  RULE  INTEGRATION  SCHEME.   THE  SUBROUTINE  WILL  ONLY  HANDLE 

C*  LINEAR  ELASTIC-PERFECTLY  PLASTIC  MATERIALS. 

C* 

C*****************************************:M**:M********************* 

c* 
c 

SUBROUTINE  PLAST(I) 
C 

IMPLICIT  REAL*8  (A-H.O-Z)C 

REAL*8  MRP.MTP.NU 

INTEGER  IZ 
SINCLUDE: 'COMMON. FOR' 
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c 

c* 

c****************************************************************** 

c* 


c» 
c* 
c* 

VARIABLES 

DSED  IN  THIS  SECTION  ARE  DEFINED  AS  FOLLOWS 

DEB 

- 

TOTAL  RADIAL  STRAIN  INCREMENT 

c* 

DET 

- 

TOTAL  TANGENTIAL  STRAIN  INCREMENT 

c* 

DERE 

- 

ELASTIC  RADIAL  STRAIN  INCREMENT 

c* 

DETE 

- 

BLASTIC  TANGENTIAL  STRAIN  INCREMENT 

c» 

DEL 

- 

DEPTH  INCREMENT 

c» 

I 

- 

RADIAL  POSITION  STATION  NUMBER 

c* 

IZ 

- 

HALF-THICKNESS  STATION  NUMBER 

c* 

DS1 

- 

RADIAL  STRBSS  INCREMENT 

c* 

DS2 

- 

TANGENTIAL  STRESS  INCREMENT 

c* 

STRT 

- 

MATRIX  HOLDING  THE  TEMPORARY  STRESSES  AT  THE 

c* 

GIVEN  LOADING  INCREMENT  ( 1 =RAD I AL ; 2=TANGENTIAL ) 

c* 

SET 

- 

TEMPORARY  EQUIVALENT  STRESS  -  USED  FOR  YIELD 

c» 

CRITERIA 

c* 

SE 

- 

EQUIVALENT  STRESS  FROM  PREVIOUS  INCREMENT 

c* 

DEPT 

- 

ARRAY  CONTAINING  THE  TEMPORARY  STRAIN  INCREMENTS 

c* 

AT  THE  GIVEN  LOADING  INCREMENT  (1=RADIAL;2= 

c* 

TANGENTIAL 

c* 

R 

- 

CONSTANT  REFERRING  TO  THE  FRACTION  BEYOND  YIBLDING 

c* 

THAT  OCCURRED  DURING  THE  COURSE  OF  THE  LAST  LOAD 

c» 

INCREMENT 

c* 

Z 

- 

DEPTH 

c* 

c* 

c 

STT(I,9] 

i=0 

.0 

STT(I,  10 )  =< 

D.O 

SSUMR=0. 

0 

SSUMT=0. 

0 

DEL=GEOf 
Z  =  0 
E=GEOM(( 

l(5)/2.0/(ADJ(6)-1.0) 

>)/GEOM(7) 

NU=GEOM(3) 

DO  500  IZ=NINT(ADJ(6) ) , 1,-1 
C 

C» 
Ctt***tt*t*t*tt*tttttt***t**t*t*t*t****t*tt*t**t*t*tt**tt****t*ttttt* 

c* 

C*     THIS  SECTION  CALCULATES  THE  STRAIN  AND  STRESS  INCREMENTS 

C*  CORRESPONDING  TO  THE  GIVEN  LOAD  INCREMENT.    IT  THEN  ADDS  THE 

C*  STRESS  INCREMENT  TO  THE  PREVIOUS  TOTAL  STRESS  TO  CALCULATE  THE 

C*  EQUIVALENT  STRESS.    THE  EQUIVALENT  STRESS  IS  COMPARED  TO  THE 

C*  UNIAXIAL  YIELD  STRESS  TO  SEE  IF  YIELDING  HAS  OCCURRED.    IF  SO, 

C*  THE  PLASTIC  STRAIN  INCREMENT  IS  SET  EQUAL  TO  THE  TOTAL  STRAIN 

C*  INCREMENT 

C* 

C*************************»****************************************** 

C* 
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DER=-STT(I,7)*(IZ-1)*DBL 
DBT=-STT(I,3)*(IZ-1)*DEL/ST(I, 1) 
DERB  =  DER-DEPTd,IZ, 1) 
DETB=DET-DEPT(I, 11,2) 
DS1=B/ ( 1 . 0-NU**2 . 0 ) * ( DERE+NU*DETE ) 
DS2=E/(1.0-NU**2.0)*(DETB+NU*DERE) 
STHT(I,  IZ,  1)=STR(I,  IZ,  D+DS1 
STRT(I, IZ,2)=STR(I,  IZ,2)+DS2 

SET(I, IZ)=SQRT(STRT(I, IZ, 1 ) **2 . O-STRT ( I , IZ, 1)*STRT(I, IZ,2) 
+  +STRT(I, IZ,2)**2) 

C 

C* 

C*      SEE  IF  YIELD  CRITERIA  IS  EXCEEDED.   IF  SO,  REZERO  THE  PLASTIC 
C*   STRAIN  INCREMENTS  AND  THE  STRBSS  INCREMENTS.   ALSO,  RETURN  TO  THE 
C*   MAIN  PROGRAM  IF  THE  OUTER  FIBER  HAS  NOT  BEGUN  TO  YIELD 
C* 
C 

IF  (SET(I.IZ)  .LT.  GE0M(6))  THEN 
DBPT(I, IZ, 1)=0.0 
DEPT(I,IZ,2)=0.0 
STT(I,9)=0.0 
STT(I, 10)=0.0 
IF  (IZ  .BQ.  NINT(ADJ(6) ) )  THEN 

RBTURN 
ENDIF 
ELSE 

SET(I,IZ)=GE0M(6) 

DEPT(I, IZ, 1)=DBR 

DEPT(I, IZ,2)=DET 

IF  (SE(I.IZ)  .LT.  GEOM(6))  THEN 

R=(SET(I, IZ)-GEOM(6) )/(SET(I, IZ)-SE(I,IZ)) 
DEPT(I,IZ,l)=R*DEPTd,IZ,l) 
DEPT(I, IZ,2)=R*DEPT(I, IZ.2) 
STRTd,  IZ,  1)=STR(I,  IZ,  1)+DS1*(1-R) 
STRT(I,IZ,2)=STHd,  IZ  ,  2  ) +DS2*  (  1-R) 
ENDIF 
ENDIF 
C 

Z=Z+DEL 
C 

500   CONTINUE 
C 
C* 

c* 

C*       PERFORM  A  TRAPAZOIDAL  RULE  INTEGRATION  OF  THE  STRAIN  OVER  THE 

C*   PLATE  HALF  THICKNESS,  AND  MULTIPLY  BY  2  TO  GIVE  THE  TOTAL  RADIAL 

C*   AND  TANGENTIAL  PLASTIC  MOMENT  INCREMENTS  (STTd.9)  AND  STTd.lO).) 

C* 

C********»:M******M*************  ****************************  ********* 

C» 

c 

Z^-DEL/2.0 

DO  540  IZ=1,NINT(ADJ(6}-1) 
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c 


Z=Z+DBL 

DMPH=(  (DKPTd,  IZ,  1)+DBPT(I,IZ+1,  1)  ) +GEOM( 3 ) * ( DBPT ( I ,  IZ.2) 
+  +DEPT(I, IZ+1,2) ) )*Z/2.0 

DMPT=((DBPT(I,IZ,2)+DBPT(I,IZ+l,2))+GKOM;3)*(DEPT(I,IZ,l) 
+  +DBPKI,  IZ+1,  1)  )  )*Z/2.0 

SSUMR=SSUMR+DMPR 

SSUMT=SSUMT+DMPT 
540   CONTINUE 

STT( I,9)=2.0*GEOM(6)/GEOM(7)/(1-GEOM(3)**2.0)*DEL*SSUMR 
STT(I, 10)=2.0*GEOM(6)/GBOM(7)/(1-GEOM(3)**2.0)*DBL*SSUMT 

WRITE(*,600)  I 
600   FORMAT(2X, ' LAT  POS  =  ',12) 

WRITE(*,700)  STT(I,9) ,STT(I,  10) 
700   FORMAT(2X, 'MRP  =  '.F20.12,'    MTP  =  '.F20.12) 

RETURN 
END 


C 

C* 

C»***************»*»»*»*********************************************»* 

c* 

C*  SUBROUTINE  UPDATE 

C* 

C*  THIS  SUBROUTINE  TAKES  THE  INCREMENTAL  STRESSES,  STRAINS, 

C*  DEFLECTIONS,  MOMENTS,  ETC  AND  ADDS  THEM  TO  THE  PREVIOUS  TOTALS 

C«  ONCE  CONVERGENCE  OF  PLASTIC  MOMENTS  HAS  OCCURRED. 

C» 

C* 

c 

SUBROUTINE  UPDATE (STT1, IMAX , OUT , LDMAT1 , LL , LDMATO ) 
C 

IMPLICIT  REAL*8  (A-H.O-Z) 

REAL*8  STT1(50,2) , LDMATO (50) ,LDMAT1(50) ,NU,DEL,E 

INTEGER  NN.OUT.LL 
SINCLUDE: 'COMMON. FOR* 
C 

NN=NINT(ADJ(2) )+l 

E=GEOM(6)/GEOM(7) 

DEL=GEOM(5)/2.0/(ADJ(6)-1.0) 

NU=GEOM' 3 ) 
C 

DO  1000  1=1, NN 

DO  905  J  =  LL,  10 

ST(I,J)=ST(I,J)+STT(I, J) 
905      CONTINUE 
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LDMATO ( I )  =  LDMAT1 ( I ) +LDMATO ( I ) 
C 

DO  980  IZ=1,NINT(ADJ(6) ) 
C 

C* 

C»     IF  THE  YIELD  CRITERIA  HAS  BEEN  EXCEBDED,  THEN  ADD  IN  THE 
C*   PLASTIC  STRAIN  INCREMENTS  AND  UPDATE  THE  TOTAL  STRESS  TO 
C*   EQUAL  THE  TEMPORARY  STRESS  CALCULATED  IN  THE  SUBROUTINE 
C*   PLAST 
C» 
C 

IF  (SE(I,IZ)  .GT.  GEOM(6))  THEN 

BP(I,IZ,1)=EP(I,IZ,1)+DEPT(I,IZ,1) 
BP(I,IZ,2)=EP(I,  IZ,2)-rDEPT(I,  IZ,2) 
STR( I, IZ, 1)=STRT(I, IZ, 1) 
STR(I, IZ,2)=STRT(I, IZ.2) 
Z=(GBOM(5)/2.0/(ADJ(6)-1.0))*(IZ-1) 
ELSE 
C 

C* 

C*     IF  THE  YIELD  CRITERIA  HAS  NOT  BEEN  EXCEEDED,  THEN  UPDATE  THE 
C*   STRESSES  AND  EQUIVALENT  STRESSES  USING  THE  ELASTIC  EQUATIONS 
C* 
C 

STR(I, IZ, 1)=-E/(1.0-NU**2.0)*(ST(I,7)*(IZ-1)*DEL- 
+  EP(I, IZ, 1)+NU*(ST(I,3)*(IZ-1)*DEL/ST(I, 1) )) 

STR(I, IZ,2)=-E/(1.0-NU**2.0)*(ST(I,3)*(IZ-1)*DEL/ST(I, 1)- 
+  BP( I, IZ,2)+NU*(ST(I,7)*(IZ-1)*DEL) ) 

SB(I, IZ)=SQRT(STR(I, IZ,1)**2.0-STR(I,IZ,1)» 
+  STR( I, IZ,2)+STR( I, IZ,2)**2) 

ENDIF 

IF  (ABS(EP(I, IZ, 1) )  .GT.  GE0M(9)  .OR.  ABS ( EP ( I , IZ , 2 ) ) 
+  .GT.  GEOM(9))  THEN 

OUT=l 
ENDIF 
980       CONTINUE 
1000  CONTINUE 
C 

LL=2 

RBTURN 

END 
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c 
c* 

Ciuuut*mnumuu*u»mnu»jnmt»*tu*tmu«iti**»*u»Uf 
C* 

C*    SUBROUTINE  BBSUL 

C* 

C*      STOBE  THE  BESULTS  AT  THE  END  OF  EVERY  10  LOAD  INCBBMENTS  IN 

C*    THE  FILE  OUT. DAT.   ALSO  SCREEN  DUMP  THESE  BESULTS. 

C* 

C**********^***********************************4^*4^***** ************ 

c* 

c 

SUBBOUTINE  BESUL(MM, STT1 , COUNT, LDMATO ) 
C 

IMPLICIT  BEAL*8  (A-H.O-Z) 
BEAL*8  STT1(50,2) ,LDMAT0(50) 
INTEGER  NN 
$INCLUDE: 'COMMON. FOR" 
C 

NN=NINT(ADJ(2))+1 
DEL=(GEOM(l)-GEOM(2))/ADJ(2) 
C 

IF  (COUNT  .EQ.  1.0)  THEN 
WRITE(15,100) 
WRITB(*,100) 
100      FORMAT(///, '************************************************' , 
+      • ********************************* ,//) 
WRITE(15, 120)  MM,LDMAT0(1)/DEL 
WRITE(*,120)  MM, LDMATO ( 1 ) /DEL 
120      FORMAT(2X, 'LOAD  INCR  =  ',13,'     LOAD  =  *,F10.4,/) 


C 


WHITE (15, 130) 
WRITE(*, 130) 
130      FORMAT(2X, 'RADPOS' ,2X, 'DEFLECT  *,2X, 'SLOPE   ',2X,'RAD  MOM  ' 
+       ,2X,'  MRP   ',2X,'   TAN  MOM   ',2X,'  MTP   ','   SHEAR    ',/) 
DO  160  1=1, NN 

STM1=ST(I,4)+ST(I,9) 
STM2=ST(I,6)+ST(I,10) 

WRITE(15,140)  ST(I,1),-ST(I,2),ST(I,3),STM1,ST(I,9), 
+  STM2,ST(I, 10) ,ST(I,5) 

WRITE (*, 140)  ST (1,1), -ST (I, 2), ST (I, 3) , STM1 , ST( I , 9) , 
+  STM2,ST(I,10),ST(I,5) 

140  FORMAT(2X,F5.2,2X,F7.5,2X,F7.4,2X,F10.0,2X,F5.0,2X,F10.0, 

+         2X.F5.0.2X.F10.0) 
160      CONTINUE 

WRITE (15, 163) 
WRITE(*,163) 
163      FOBMAT(/,2X, ' 

+      ' ') 

DO  1000  1=1, NN 

IF  (ST(I,9)  .NE.  0.0  .08.  ST(I,10)  .NE.  0.0)  THEN 
WHITE (15, 170)  ST(I.l) 
WBITE(*,170)  ST(I.l) 
170  FOBMAT(/,2X, 'BADIAL  POSITION  =  *,F5.2) 
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WBITE(15,165) 
WRITE(*,165) 
165  F0HMAT(/,2X, '  DEPTH   ',2X,'    BBP      ',2X,'   SHAD 

+  2X,*    ETP     ',2X,'   STAN    ',2X,'    SB     ',/) 

DO  200  IZ=1,NINT(ADJ(6) ) 

IF  (SB(I,IZ)  .GT.  GEOM(6))  THEN 

Z=(GBOM(5)/2.0/(ADJ(6)-1.0))*(IZ-1) 

WRITB(15,180)  Z,EP(I,IZ,1),STR(I,IZ,1),BP(I,IZ,2), 
+  STR(I,IZ,2),SE(I,IZ) 

WRITE (*, 180)  Z,EP(I,IZ,1),STR(I,IZ,1),BP(I,IZ,2), 
+  STR(I,IZ,2),SE(I,IZ) 

180  FOHMAT(2X,F6.4,3X,F10.8,3X,F8. 1 , 3X , F10 . 8 , 3X , F8. 1, 

+  3X.F8.1) 

BNDIF 
200  CONTINUE 

ENDIF 
1000     CONTINUE 
ENDIF 
RETURN 
END 


COMMON  GE0M(9) ,ADJ(7) ,ST(50, 10) 

COMMON  /PLSCM/  EP(50,50,2),STT(50,10),DEPT(50,50,2),STR(50,50,2), 
+     SB (50, 50) ,STRT(50,50,2) , SET (50, 50) 
REAL*8  ADJ,ST,EP,STT,DEPT,STR,SE,STRT,SET 


132 


LIST  OF  REFERENCES 


1.   Stern,  M.,  "A  General  Boundary  Integral  Formulation  for 
Plate  Bending  Problems",  The  International  Journal  of 
Solids  Structures,  15,  pp.  769-782  (1979) 


2.  Bezine,  G.,  "Boundary  Integral  Formulation  for  Plate 
Flexure  With  Arbitrary  Boundary  Conditions",  Mechanics 
Research  Communications,  5,  pp.  197-206  (1978) 

3.  Moshaiov,  A.  and  Vorus,  W.  S.,  "Elasto-Plast ic  Plate 
Bending  Analysis  by  a  Boundary  Element  Method  With  Initial 
Plastic  Moments",  The  International  Journal  of  Solids  and 
Structures ,  article  in  print  (1986) 

4.  Karaiya,  N.  and  Sawaki,  Y.,  "An  Integral  Equation  Approach 
to  the  Finite  Deflection  of  Elastic  Plates",  Journal  of 
Non-Linear  Mechanics,  17,  pp. 187-194  (1982) 


133 


5.  Butterfield,  R. ,  "New  Concepts  Illustrated  by  Old 
Problems",  Developments  in  Boundary  Element  Methods, 
ed.  P  Banerjee  and  R.  Butterfield,  Applied  Science 
Publishers  (1980) 

6.  Timoshenko,  S.  A.  and  Woinowsky-Kreiger ,  S.,  Theory  of 
Plates  and  Shells,  2nd  ed.,  McGraw  Hill,  New  York  (1959) 

7.  Tuma,  Jan  J. ,  Theory  and  Problems  of  Structural  Analysis, 
McGraw  Hill,  New  York  (1969) 

8.  Mendelson,  A.,  Plasticity:  Theory  and  Application, 
MacMillan  Co.,  New  York  (1968) 

9.  Roark,  R.  and  Young,  W.,  Formulas  for  Stress  and  Strain, 
5th  ed. ,  MacGraw  Hill,  New  York  (1975) 

10.  Armen,  H.  Jr.,  Pifko,  A.  and  Levine,  H.,  "Finite  Element 
Analysis  of  Structures  in  the  Plastic  Range",  NASA  Report 
CR-1649,  National  Aeronautics  and  Space  Administration, 
Washington,  DC  (1971) 


134 


DUDLE 

r^f;  'HOOL 

m^  BBIA  93945  300* 


Eareckson 

Elasto-plastic  analy- 

loaded  circular  plates 
"sxng  Green's  functions. 


Thesis 

E1146 

c.l 


2!92<t9 


Eareckson 


Elasto-plastic  analy- 
sis of  axisymmetrically 
loaded  circular  plates 
using  Green's  functions. 


